■  jju  ■  Defence  Research  and  Recherche  et  developpement 

■  ■  Developmenf  Canada  pour  la  defense  Canada 


Advanced  beamformers  for  3D  ultrasound 
systems  deploying  linear  and  planar 
phased  array  probes 


Stergios  Stergiopoulos 
Amar  Dhanantwari 


DISTRIBUTION  STATEMENT  A 

Approved  for  Public  Release 
Distribution  Unlimited 


Defence  R&D  Canada  -  Toronto 

Technical  Report 
DRDC  Toronto  TR  2002-058 
May  2002 


20030129  m 


Canada 


Advanced  beamformers  for  3D 
ultrasound  systems  deploying  linear  and 
planar  phased  array  probes 


Stergios  Stergiopoulos 
Amar  Dhanantwari 


Defence  R&D  Canada  -  Toronto 

Technical  Report 
DRDC  Toronto  TR  2002-058 
May  2002 


Stergi 


deigiopoulos,  Ph.D. 


Amar  Dhanantwari,  Ph.D. 


Approved  for  release  by 

1  — . 

K.M.  Sutton 

Chair,  Document  Review  and  Library  Committee 


Her  Majesty  the  Queen  as  represented  by  the  Minister  of  National  Defence,  2002 
Sa  majeste  la  reine,  representee  par  le  ministre  de  la  Defense  nationale,  2002 


Abstract 


The  present  report  outlines  the  development  of  an  advanced  digital  ultrasound  imaging 
technology  leading  to  a  next-generation  field-portable  4D  diagnostic  imaging  system.  The 
proposed  technology  consists  of: 

•  Adaptive  beamforming  algorithms  to  provide  high  image-resolution  for 
ultrasound  diagnostic  systems; 

•  3D  visualization  imaging  techniques  to  assist  imaging  during  field-deployable 
minimally  invasive  operations; 

•  time-reversal  pulse  design  for  ultrasound  signals  to  eliminate  aberration 
effects,  which  cause  fuzziness  in  reconstructed  ultrasound  images;  and 

•  investigations  on  computing  architectures  and  planar  ultrasound  sensor  arrays 
allowing  system  integration  of  the  proposed  technologies  as  a  compact- 
portable  system  for  field-deployable  operations. 

The  aim  is  to  investigate  a  new  ultrasonic  sensing  and  imaging  technology  that  would  lead 
to  the  design  of  portable,  compact  and  field  deployable  3D  ultrasound  diagnostic  systems. 
The  performance  characteristics  of  the  proposed  diagnostic  systems  include  minimization  of 
the  aberration  effects  and  significant  improvement  in  image  resolution  to  allow  for  tissue 
identification  and  3D  tomographic  imaging  of  internal  body-organs.  The  adaptive  processing 
schemes  of  this  study  are  characterized  as  Dual-Use  technologies  and  have  been  tested  in 
operational  active  sonars  of  the  Canadian  Navy.  Real  data  results  have  shown  that  they 
provide  array  gain  improvements  for  signals  embedded  in  anisotropic  noise  fields,  similar  to 
those  in  the  human  body. 
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Resume 


Le  present  rapport  traite  de  la  mise  au  point  d’une  technologie  d’imagerie  a  ultrasons 
numerique  evoluee  en  vue  de  creer  un  systeme  d’imagerie  4D  de  prochaine  generation 
pouvant  etre  utilise  sur  place  a  des  fins  de  diagnostic.  La  technologie  proposee  comprend  les 
travaux  suivants  : 

•  techniques  de  formation  de  faisceaux  adaptatives  [1-6]  permettant  de  creer  des  images 
a  haute  resolution  pour  les  systemes  de  diagnostic  a  ultrasons; 

•  techniques  d’imagerie  3D  afin  de  faciliter  l’imagerie  durant  les  operations  avec 
effraction  minimale  effectuees  sur  place; 

•  concept  d’impulsions  a  inversion  temporelle  pour  les  signaux  ultrasons  afin 
d’eliminer  les  effets  d’aberration,  qui  rendent  floues  les  images  reconstruites  des 
echographies; 

•  etudes  sur  les  architectures  informatiques  et  les  reseaux  plans  de  capteurs  d’ultrasons 
afin  de  permettre  l’integration  des  technologies  proposees  dans  un  systeme  portatif 
compact  utilisable  pour  les  operations  sur  place. 

Notre  objectif  consiste  a  etudier  une  nouvelle  technologie  de  detection  et  d'imagerie 
ultrasoniques  qui  permettrait  de  concevoir  des  systemes  a  ultrasons  3D  portatifs,  compacts  et 
utilisables  sur  place  aux  fins  de  diagnostic.  Les  caracteristiques  de  performance  des 
systemes  de  diagnostic  proposes  comprendraient  la  reduction  des  effets  d’aberration  et  une 
amelioration  considerable  de  la  resolution  des  images  afin  de  permettre  1 ’identification  des 
tissus  et  l’imagerie  tomographique  3D  des  organes  internes  du  corps.  Les  techniques  de 
traitement  adaptatif  utilisees  dans  la  presente  etude  sont  appelees  technologies  a  double 
usage  et  ont  ete  mises  a  l’essai  dans  des  sonars  actifs  operationnels  de  la  Marine  canadienne 
[1,4-6],  Les  resultats  obtenus  en  situation  reelle  ont  demontre  que  ces  techniques  permettent 
d’ameliorer  dans  les  reseaux  le  gain  des  signaux  incorpores  dans  des  champs  sonores 
anisotropes  similaires  a  ceux  qu’on  retrouve  dans  le  corps  humain. 
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Executive  summary 


The  aim  of  this  report  is  to  bring  together  some  of  the  most  recent  theoretical  developments 
on  ultrasound  beamformers;  and  provide  suggestions  of  how  modem  technology  can  be 
applied  to  the  development  of  current  and  next  generation  3D  ultrasound  systems.  It  will 
focus  on  the  development  of  an  advanced  beamforming  structure  that  allows  the 
implementation  of  adaptive  and  synthetic  aperture  signal  processing  techniques  in 
ultrasound  systems  deploying  multi-dimensional  arrays  of  sensors. 

The  primary  objective  is  to  develop  an  Ultrasonic  Sensing  and  Imaging  Technology  that  will 
facilitate  non-invasive  medical  imaging  for  treatment  assessment  of  injured  CF  personnel, 
combat  casualties,  and  victims  rescued  during  search  and  rescue  operations.  The 
requirements  for  non-invasive  medical  diagnostic  procedures  can  be  addressed  with  the 
development  of  compact,  portable  and  field-deployable  3D  ultrasound  medical  imaging 
systems.  Because  imaging  is  of  critical  importance  in  facilitating  rapid  diagnosis  and 
therapy,  the  major  manufacturers  of  imaging  equipment  are  introducing  CT  &  MRI 
machines  that  allow  surgical  procedures  to  be  performed  under  direct  imaging  using 
minimally  invasive  techniques.  However,  these  imaging  approaches  are  expensive  and 
would  occupy  a  complete  surgical  suite.  Unlike  CT  or  MRI,  which  are  often  located  in 
dedicated  suites  specifically  built  to  accommodate  their  needs,  ultrasound  is  a  much  cheaper 
and  more  flexible  imaging  modality  that  can  be  deployed  in  remote  and  harsh  operational 
areas  of  interest  to  the  CF.  Moreover,  in  contrast  with  clinical  CT  and  MRI,  ultrasound 
systems  allow  the  surgeon  or  therapist  to  see  what  is  happening  as  it  occurs,  in  real  time. 
However,  even  the  most-advanced  state-of-the-art  medical  ultrasound  imaging  systems 
suffer  from  very  poor  image  resolution,  which  is  the  result  of  the  very  small  size  of  deployed 
arrays  of  sensors  and  the  distortion  effects  caused  by  the  human-body’s  non-linear 
propagation  characteristics.  In  summary,  the  existing  ultrasound  imaging  technology  is  not 
field-deployable  and  incapable  of  delivering  high-quality  images. 

The  development  of  an  advanced  4-Dimensional  (4D,  i.e.  3D  +  Time)  ultra-sound  system 
technology,  detailed  in  the  present  report,  addresses  ultrasound  requirements  for  high-image 
resolution  capabilities  with  small  size  arrays  for  medical  diagnostic  and  non-destructive 
system  applications.  The  present  investigation  has  been  supported  by  the  Technology 
Investment  Fund  of  Defence  R&D  Canada  (DRDC)  and  two  European-Canadian 
collaborative  projects  entitled:  “MITTUG:  (Minimally  Invasive  Tumour  Therapy  3D 
Ultrasound  Guided)”  and  “ADUMS:  (Adaptive  4D  Processing  for  Digital  Ultrasound 
Medical  Systems)”  that  have  been  awarded  with  funds  from  the  European  Commission  Fifth 
Framework  1ST  program. 

In  summary,  it  is  anticipated  that  the  results  of  our  4D  ultra-sound  system  technology  would 
provide  a  major  support  in  terms  of  innovation  to  the  emerging  Canadian  ultra-sound 
industrial  sector,  since  the  anticipated  high  image  resolution  capabilities  can  be  exploited  by 
other  medical  imaging  diagnostic  and  non-destructive  system  applications. 
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Sommaire 


Le  present  rapport  vise  a  rassembler  certaines  des  avancees  theoriques  les  plus  recentes  dans 
le  domaine  des  conformateurs  de  faisceaux  ultrasonores  et  de  suggerer  des  faqons  dont  la 
technologie  modeme  peut  etre  appliquee  au  developpement  d’appareils  a  ultrasons  3D  de  la 
generation  actuelle  et  de  la  prochaine  generation.  Le  rapport  portera  principalement  sur  la 
mise  au  point  d’une  structure  evoluee  de  conformation  de  faisceaux  afin  de  permettre 
l’utilisation  de  techniques  de  traitement  d’ouverture  synthetique  et  de  traitement  adaptatif  de 
signal  dans  des  appareils  a  ultrasons  deployant  des  reseaux  multidimensionnels  de  capteurs. 

Le  principal  objectif  consiste  a  developper  une  technologie  de  detection  et  d’imagerie 
ultrasoniques  qui  facilitera  l’application  de  l’imagerie  medicale  non  effractive  a  revaluation 
medicale  et  au  traitement  du  personnel  blesse  des  FC,  des  blesses  au  combat  et  des  victimes 
secourues  durant  les  operations  de  recherche  et  de  sauvetage.  Les  besoins  en  matiere  de 
procedures  de  diagnostic  medical  non  effractives  peuvent  etre  satisfaits  grace  a  la  mise  au 
point  de  systemes  d’imagerie  medicale  a  ultrasons  3D  compacts,  portatifs  et  utilisables  sur 
place.  Etant  donne  que  l’imagerie  revet  une  importance  critique  dans  l’etablissement  rapide 
d’un  diagnostic  et  d’un  traitement,  les  principaux  fabricants  d’equipement  d’imagerie  ont 
introduit  des  appareils  de  tomodensitometrie  et  d’imagerie  par  resonance  magnetique  (IRM) 
qui  permettent  de  realiser  des  interventions  chirurgicales  a  l’aide  de  techniques  d’imagerie 
directe  avec  effraction  minimale.  Cependant,  ces  techniques  d’imagerie  sont  dispendieuses  et 
necessiteraient  l’utilisation  d’une  salle  de  chirurgie  complete.  Contrairement  aux  appareils  de 
tomodensitometrie  et  d’IRM,  qui  sont  souvent  installes  dans  des  salles  specialises  conques 
specialement  en  fonction  de  leurs  besoins,  l’echographie  est  une  methode  d’imagerie  bien 
moins  chere  et  plus  souple  qui  peut  etre  deployee  dans  des  secteurs  operationnels  eloignes  et 
hostiles  qui  interessent  les  FC.  De  plus,  a  1’ inverse  des  tomodensitometres  ou  des  appareils 
d’IRM  cliniques,  les  systemes  a  ultrasons  permettent  au  chirurgien  ou  au  therapeute  de  voir 
immediatement  ce  qui  se  produit,  en  temps  reel.  Toutefois,  meme  dans  les  systemes 
d’imagerie  medicale  a  ultrasons  evolues,  a  la  fine  pointe,  la  resolution  des  images  est  faible  en 
raison  de  la  tres  petite  taille  des  reseaux  de  capteurs  deployes  et  des  effets  de  distorsion  causes 
par  les  caracteristiques  de  propagation  non  lineaires  du  corps  humain.  En  resume,  la 
technologie  d’imagerie  a  ultrasons  existante  ne  peut  pas  etre  deployee  sur  place  et  ne  permet 
pas  de  foumir  des  images  de  grande  qualite. 

Le  developpement  d’une  technologie  d’echographie  4D  (3D  +  temps)  evoluee,  presentee  en 
detail  dans  le  present  rapport,  repond  aux  besoins  en  matiere  de  haute  capacite  de  resolution 
des  images  et  de  petite  taille  des  reseaux  pour  des  utilisations  dans  le  domaine  du  diagnostic 
medical  et  dans  des  systemes  non  destructeurs.  La  presente  etude  a  requ  l’appui  du  Fonds 
d’investissement  technologique  de  RDDC  et  de  deux  projets  intitules  :  MITTUG  (traitement 
de  tumeurs  avec  effraction  minimale  sous  echoguidage  3D)  et  ADUMS  (traitement  adaptatif 
4D  pour  systeme  d’echographie  medicale  numerique)  issus  de  la  collaboration  entre  le  Canada 
et  l’Europe  et  finances  par  le  Cinquieme  programme  cadre  1ST  de  la  Commission  europeenne. 

En  resume,  on  prevoit  que  les  resultats  de  notre  etude  sur  la  technologie  d’echographie  4D 
serviront  de  fondement  important  a  l’innovation  dans  le  secteur  industriel  canadien  des 
dispositifs  a  ultrasons  en  emergence,  etant  donne  que  les  capacites  d’imagerie  a  haute 
resolution  prevues  pourront  etre  exploitees  dans  d’autres  applications  d'imagcrie  destinees 
au  diagnostic  medical  et  a  des  systemes  non  destructeurs. 
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1.  Background 


In  general,  the  mainstream  conventional  signal  processing  of  ultrasound  systems  consists  of 
a  selection  of  temporal  and  spatial  processing  algorithms  [1,60].  These  algorithms  are 
designed  to  increase  the  signal-to-noise  ratio  for  improved  signal  delectability  while 
simultaneously  providing  parameter  estimates  such  as  Doppler  and  bearing  for  incorporation 
into  the  image  reconstruction  process.  Their  implementation  in  real  time  systems  had  been 
directed  at  providing  high-quality,  artifact-free  conventional  beamformers,  currently  used  in 
operational  ultrasound.  However,  aberration  effects  associated  with  ultrasound  system 
applications  suggest  that  fundamentally  new  concepts  need  to  be  introduced  into  the  signal 
processing  structure  of  next-generation  ultrasound. 

To  provide  a  context  for  the  material  contained  in  this  report,  it  would  seem  appropriate  to 
review  briefly  the  beamforming  functionality  of  a  conventional  3D  ultrasound  system  and  to 
contrast  this  functionality  with  the  basic  requirements  of  high-performance  3D  systems 
deploying  multi-dimensional  arrays  of  sensors. 


1.1  Current  State-of-the-Art  of  3D  Ultrasound  System 
Technology 

The  requirements  for  non-invasive  and  minimally  invasive  medical  diagnostic  procedures 
can  be  addressed  with  the  development  of  high-resolution  3D  ultrasound  medical  imaging 
systems.  However,  even  the  most-advanced  state-of-the-art  medical  ultrasound  imaging 
systems  suffer  from  very  poor  image  resolution,  which  is  the  result  of  the  very  small  size  of 
deployed  arrays  of  sensors  and  the  distortion  effects  caused  by  the  human-body’s  non-linear 
propagation  characteristics.  In  particular,  some  of  the  limitations  (e.g.  resolution)  of 
ultrasound  imaging  are  related  to  fundamental  physical  aspects  of  the  ultrasound  transducer 
and  the  interaction  of  ultrasound  with  tissues  (e.g.  aberration  effects).  In  addition  to 
fundamental  limitations  are  limitations  related  to  the  display  of  the  ultrasound  image  in  an 
efficient  manner  allowing  the  physician  to  extract  relevant  information  accurately  and 
reproducibly.  Specifically,  in  this  report  we  will  address  the  following  limitations: 

•  Conventional  ultrasound  images  are  2D,  hence,  the  physician  must  mentally 
integrate  multiple  images  to  develop  a  3D  impression  of  the 
anatomy/pathology  during  procedure.  This  practice  is  time-consuming, 
inefficient,  and  requires  a  highly  skilled  operator,  all  of  which  can  potentially 
lead  to  incorrect  diagnostic  and  therapeutic  decisions. 

•  Often  the  physician  requires  accurate  estimation  of  tumour  and  organ  volume. 
The  variability  in  ultrasound  imaging  and  volume  measurements  using  a 
conventional  2D  technique  is  high,  because  current  ultrasound  volume 
measurement  techniques  assume  an  idealized  elliptical  shape  and  use  only 
simple  measures  of  the  width  in  two  views  [64].  3D  images  will  provide 
means  to  obtain  accurate  and  precise  organ  and  tumor  volume  estimates  [65], 
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•  It  is  difficult  to  localize  the  thin  2D  ultrasound  image  plane  in  the  organ,  and 
difficult  to  reproduce  a  particular  image  location  at  a  later  time,  making  2D 
ultrasound  a  limited  imaging  modality  for  monitoring  of  disease 
progression/regression  and  follow-up  patient  studies. 

In  summary,  on  issues  of  3D  visualization,  we  propose  to  overcome  these  limitations  by 
combining  the  proposed  advanced  3D  beamfoming  structure  defined  in  the  present  report 
with  existing  3D  visualization  display  technologies  [60,  63-68],  To  contrast,  however,  the 
proposed  3D  ultra-sound  technology  development,  it  is  essential  to  review  the  current  state- 
of-the-art  of  ultrasound  technology  on  3D  visualization  and  beamforming  processing. 

1.1.1  Current  Beamforming  Structure  of  Ultrasound  Systems 

A  state-of-the-art  transducer  array  of  an  operational  ultra-sound  system  is 
either  linear  or  curvilinear  depending  on  the  application;  and  for  each 
deployed  array,  the  number  of  transducers  is  in  the  range  of  96  to  256 
elements.  However,  only  a  small  number  of  transducers  of  a  given  array  are 
beamformed  coherently  to  reconstruct  a  tomography  image  of  interest. 
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Figure  1.  Typical  beamformers  for  linear  arrays  of  ultrasound  systems.  The  shaded  sub-apertures 
indicate  the  selected  spatial  locations  for  the  formation  of  synthetic  aperture  according  to  ETAM 

algorithm. 


A  typical  number  of  transducers  that  may  be  included  in  the  beamforming 
structure  of  an  ultrasound  system  is  in  the  range  of  32  to  128  elements.  Thus, 
the  system  array  gain  may  be  smaller  by  approximately  10xLog10(3)  ~  5  dBs 
than  that  available  by  the  deployed  array.  Figure  1  illustrates  the  basic 
processing  steps  associated  with  the  ultrasound  beamforming  process.  For  the 
sake  of  simplicity  and  without  any  loss  of  generality,  the  array  in  Figure  1  is 
considered  to  be  linear  with  96-elements  that  can  be  used  to  transmit  and 
receive  the  radiated  ultrasound  field.  As  shown  in  Figure  1,  for  each  active 
transmission,  the  ultrasound  beamformer  processes  coherently  the  received 
signal  of  32-elements  only,  which  is  a  sub-aperture  of  the  96-element 
deployed  array.  The  active  transmission  takes  place  approximately  every  t  = 
0.3ms  ,  depending  on  the  desired  penetration  depth  in  the  body.  The  beam 
steering  process  is  at  the  broadside. 


When  an  active  transmission  is  completed,  the  receiving  32-element  sub¬ 
aperture  is  shifted  to  the  right,  as  shown  in  Figure  1.  Thus,  to  make  use  of  all 
the  96-elements  of  the  deployed  array,  the  32-element  beamforming  process 
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is  repeated  64  times,  generating  64  broadside  beams.  In  other  words,  it  takes 
approximately  64x0.3  ms  =  20  ms  to  reconstruct  a  2-D  tomography  image  of 
interest.  As  a  result,  the  resolution  characteristics  of  the  reconstructed  image 
are  defined  by  the  array  gain  of  the  beamformer  and  the  temporal  sampling  of 
the  beam  or  element  time  series,  for  analog  or  digital  beamformers, 
respectively.  In  the  specific  case  of  Figure  I,  the  pixel  resolution  along  the 
horizontal  x-axis  of  a  reconstructed  tomography  image  is  defined  by  the 
angular  resolution  along  azimuth  of  the  32-element  beamformer.  This 
resolution  is  usually  being  improved  by  means  of  interpolation,  which  defines 
the  basic  difference  between  beamformers  of  different  ultrasound  systems. 

The  pixel  resolution  along  the  vertical  y-axis  of  the  reconstructed  image  is 
defined  by  the  sampling  rate,  which  is  always  very  high  and  it  is  not  a  major 
concern  in  ultrasound  system  applications.  Thus,  improvements  of  image 
resolution  in  ultrasound  system  applications  requires  mainly  higher  angular 
resolution  or  very  narrow  beamwidth,  which  means  longer  arrays  and  longer 
sub-apertures  for  the  beamforming  process  with  consequent  technical  and 
operational  implications  and  higher  system  manufacturing  cost. 


Figure  2.  Tomography  images  of  electronic  chip  obtained  with  different  ultrasound  systems.  Left 
image  represents  output  of  a  low-end  ultra-sound  system  with  a  linear  array  of  96-sensor  and  32- 
sensor  beamformer  with  frequency  range  centered  at  7 -MHz  with  6-MHz  bandwidth.  Right  image 
was  provided  by  high-end  cardiac  ultrasound  system  deploying  a  curvilinear  128-sensor  array  with 
48-sensor  beamformer  with  frequency  range  centered  at  3.5-MHz  with  3-MHz  bandwidth. 

To  emphasize  the  limitations  of  the  current  ultrasound  technology,  Figure  2 
provides  two  tomography  images  of  an  electronic  chip  that  have  been 
obtained  with  two  different  ultrasound  systems.  The  image  at  the  left  of 
Figure  2  represents  the  output  of  a  low-cost  ultra-sound  system  with  a  linear 
array  of  96  elements  and  with  a  32  element  beamformer  at  a  frequency  range 
centered  at  7  MHz  with  6  MHz  bandwidth.  The  cost  of  this  system  was 
$30,000.  The  image  at  the  right  of  Figure  2  was  provided  by  a  high-end 
cardiac  ultrasound  system  deploying  a  curvilinear  128-element  array  with  a 
48-element  beamformer  at  a  frequency  range  centered  at  3.5  MHz  with  3 
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MHz  bandwidth.  The  cost  of  this  system  was  $280,000.  Since  the 
beamforming  characteristics  for  both  these  two  ultrasound  systems  are 
similar,  the  resolution  characteristics  of  the  reconstructed  images  of  Figure  2 
are  very  similar,  as  expected. 

The  main  advantages  of  this  simplified  beamforming  structure  are  the 
following: 

The  generation  of  broadside  beams  allows  the  use  of  frequency  regimes  that 
are  higher  than  the  corresponding  spatial-aliasing  frequency  of  the  sensor 
spacing  of  the  ultrasound  probe.  This  is  because  side-lobe  artifacts  due  to 
spatial  aliasing  are  insignificant  for  beams  with  broadside  beamsteering.  For 
example,  the  probe  that  was  used  to  generate  the  left  image  in  Figure  2  with 
7-MHz  center  frequency  had  0.35mm  sensor  spacing  that  corresponds  to  2.2 
MHz  design  frequency  for  phased  array  processing  to  avoid  spatial  aliasing 
artifacts. 

The  advantage  (suppression  of  spatial-aliasing  artifacts)  provided  by  the 
broadside  beam-steering  process  has  been  used  effectively  by  illumination 
techniques  using  higher  order-harmonics  to  achieve  deeper  penetration  with 
corresponding  higher  image  resolution  along  the  temporal  axis. 

The  field  of  view  of  the  probe  may  be  larger  than  the  aperture  size  required 
by  the  broadside  beamformer.  This  approach  minimizes  the  hardware 
complexity  of  the  analog-to-digital  converter  (A/DC)  by  using  a  multiplexer 
to  control  the  data  acquisition  process  of  a  probe  with  a  larger  number  of 
sensors  than  those  being  used  by  the  broadside  focus  beamformer. 

Until  recently,  the  above  innovative  approaches  have  served  well  the 
ultrasound  system  requirements  by  providing  practical  alternatives  to 
technical  problems  that  were  due  mainly  to  limitations  in  the  maximum 
number  of  channels  deployed  by  A/DC  units  and  the  limited  capabilities  of 
computing  architectures.  Presently,  these  type  of  technology  limitations  do 
not  exist.  Thus,  new  technology  development  options  have  become  available 
to  exploit  the  vast  experience  from  phased  array  beamformers  that  have  been 
advanced  by  the  sonar  and  radar  research  communities. 

In  fact,  the  introduction  of  linear  phased  array  probes  for  cardiac  applications 
is  the  first  successful  step  toward  this  direction. 

The  use,  however,  of  linear  arrays  introduces  another  major  problem  in  terms 
of  false  targets,  a  problem  that  has  been  identified  by  both  the  ultrasound  and 
sonar  researchers  using  towed  sonar  arrays  [60].  In  particular,  a  linear  array 
provides  angular  resolution  within  the  tomography  plane  (B-scan)  that  the 
beam  steering  is  formed.  The  angular  resolution,  however,  of  the  beam 
steering  vectors  of  linear  arrays  is  omnidirectional  in  the  plane  perpendicular 
to  the  B-scan  plane.  Thus,  reflections  from  surrounding  organs  cannot  be 
spatially  resolved  by  the  steered  beams  of  a  line  array;  and  they  appear  as 
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false  targets  in  towed  array  sonars,  or  false  components  of  a  reconstructed 
image  by  a  linear  ultrasound  probe. 

To  address  the  problem  of  false  components  in  the  reconstructed  image,  the 
1.5D  and  1.75D  ultrasound  array  probes  have  been  introduced  that  consist  of 
linear  arrays  stacked  as  partially  planar  arrays.  In  particular,  the  GE  1 .75D 
ultrasound  array  probe  consists  of  8  linear  arrays  with  128-sensors  each  and 
with  0.2mm  sensor  spacing.  The  linear  array  spacing  is  1 .5mm.  Moreover,  the 
1.5D  array  probes  consist  of  4  linear  arrays  with  128-sensors  each. 

Thus,  the  steered  beams  for  both  the  1 ,5D  and  1 .75D  type  of  array  probes  are 
3-dimensional  (3D)  and  have  the  property  to  resolve  the  angular  components 
of  ultrasound  reflected  signals  along  azimuth  and  elevation.  Although  the  3D 
beams  of  1.75D  arrays  may  be  viewed  as  the  first  step  for  3D  ultrasound 
imaging,  they  do  not  have  sufficient  angular  resolution  capabilities  along 
elevation  to  generate  3D  ultrasound  volumes.  The  proposed  3D  beamforming 
structure  and  the  relevant  3D  ultrasound  system  development  based  on  a 
(32x32)  planar  array  probe  attempt  to  address  the  above  limitations  and  the 
results  are  discussed  in  the  present  report.  However,  at  this  point,  it  is 
considered  appropriate  to  briefly  review  the  current  state  of  the  art  in  3D 
ultrasound  technology. 

1.1.2  Current  Technology  Concept  of  3D  Visualization 
Methods  for  Ultrasound  Systems 

Current  3D  Ultrasound  imaging  systems  have  3  components:  image 
acquisition,  reconstruction  of  the  3D  image,  and  display  [66,67],  The  first 
component  is  crucial  in  ensuring  that  optimal  image  quality  is  achieved.  In 
producing  a  3D  image,  the  conventional  transducer  is  moved  over  the 
anatomy  while  2D  images  are  digitized  and  stored  in  a  microcomputer,  as 
shown  in  Figure  3.  To  reconstruct  the  3D  geometry  without  geometric 
distortion,  the  relative  position  and  angulation  of  the  acquired  2D  images 
must  be  known  accurately.  Over  the  past  6  years  there  have  been  numerous 
developments  and  evaluating  techniques  for  obtaining  3D  ultrasound  images 
using  two  approaches:  mechanical  scanning  and  freehand  scanning  [65,66,68- 

71]- 
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Figure  3. Mechanical  scanning  of  ultrasound  probes  for  image 
acquisition  of  2D  B-scans  to  obtain  3D  ultrasound  images 
through  volume  rendering. 

Mechanical  Scanning:  Based  on  earlier  work,  Fraunhofer  (National 
Research  Center,  Darmstadt,  Germany)  and  the  Robarts  Research  Institute 
(RRI  of  the  University  of  Western  Ontario,  London  Canada),  have  developed 
systems  in  which  the  ultrasound  transducer  is  mounted  in  a  special  assembly, 
which  can  be  driven  by  a  motor  to  move  in  a  linear  fashion  over  the  skin  or 
tilted  in  equal  angular  steps.  The  movement  can  be  continuous,  or 
intermittent  at  cardiac  and/or  respiratory  intervals  [65,70].  In  addition,  the 
spatial-sampling  frequency  of  the  image  acquisition  can  be  adjusted  based  on 
the  elevational  resolution  of  the  transducer  and  the  depth  of  the  region-of- 
interest.  For  linear  scanning,  they  collect  140  images  (336  x  352  pixels  each) 
at  0.5  mm  intervals  in  a  time  that  depends  on  the  ultrasound  machine  frame 
rate  and  whether  cardiac  gating  is  used.  All  scanning  parameters  can  be 
adjusted  depending  on  the  experiment  and  type  of  acquisition.  For  example 
for  3D  B-mode  -  they  typically  use  2  or  3  focal  zones  resulting  in  about  15 
ffames/sec  and  a  total  3-D  scanning  time  of  9  sec  for  140  images. 

Freehand  Scanning:  Although  the  mechanical  scanning  approach 
produces  accurate  3D  ultrasound  images,  the  mechanical  assembly  is  bulky 
and  not  convenient  for  the  operator.  In  addition,  the  mechanical  constraint 
does  not  permit  its  use  in  imaging  larger  structures  such  as  the  liver  or  fetus. 
Thus,  alternative  techniques  have  been  introduced  that  maintain  the  flexibility 
of  the  2D  exams  yet  produce  3D  images.  These  techniques  include  freehand 
scanning  systems,  in  which  a  magnetic  positioning  and  orientation 
measurement  (POM)  device  is  mounted  on  the  transducer  [63,68,71-73],  To 
produce  a  3D  image,  the  operator  manually  moves  the  hand-held  transducer, 
while  the  POM  device  transfers  the  position  and  orientation  coordinates  of 
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the  transducer  to  a  microcomputer.  At  the  same  time,  2D  images  are  digitized 
by  the  same  computer  and  associated  with  the  appropriate  coordinates.  After 
the  necessary  number  of  2D  images  are  acquired  (typically  60  -  1 60),  the 
computer  reconstructs  the  3D  image.  Care  is  taken  to  scan  the  patient 
sufficiently  slowly  so  that  the  region  of  interest  is  scanned  with  no  gaps. 
Typically,  the  scan  lasts  4-11  sec.  while  the  patient  holds  its  breath. 
Although  this  technique  does  produce  useful  images,  it  still  suffers  form 
major  limitations  that  precludes  its  use  for  general  diagnostic  procedures. 

Most  importantly,  the  manual  scanning  of  the  3D  space  with  a  linear  array 
does  not  eliminate  the  false  components  of  the  reconstructed  B-scan  images 
that  were  discussed  in  the  previous  section. 


1.2  Next  Generation  Real-Time  3D  Ultarsound  System 
Technology 

As  part  of  our  DRDC’s  Technology  Investment  Fund  (TIF)  project  objectives,  the 
development  of  an  ultrasound  imaging  technology  leading  to  a  next-generation  high- 
resolution  real  3D  diagnostic  imaging  system  has  been  initiated.  The  proposed  technology 
includes: 

•  Synthetic  aperture  processing  to  sample  large  size  2D  arrays  for  3D  phased 
array  beamforming, 

•  adaptive  beamforming  to  effectively  increase  the  angular  resolution 

•  32x32-sensor  phased  planar  arrays  with  uniform  0.4mm  sensor  spacing  to 
achieve  maximum  sensitivity. 

•  Integration  of  Fraunhofer’s  3D  and  4D  visualization  schemes  with  the  image 
reconstruction  process  of  the  3D  ultrasound  beamformer. 

As  stated  above,  the  long  term  objective  of  the  present  development  is  to  replace  the  existing 
beamformers  of  ultrasound  systems  with  an  adaptive  beamforming  staicture  in  order  to 
maximize  the  array  gain  and  image  resolution  capabilities  of  2D  and  3D  ultrasound  systems. 

1.2.1  Synthetic  Aperture  Processing  in  the  Beamforming 
Structure  of  Ultrasound  Systems 

At  this  point  it  is  important  to  review  a  few  fundamental  physical  arguments 
associated  with  synthetic  aperture  processing  for  sonar  and  ultrasound 
systems.  In  the  past  [13]  there  was  a  conventional  wisdom  regarding 
synthetic  aperture  techniques,  which  held  that  practical  limitations  prevent 
them  from  being  applicable  to  real-world  systems.  The  issues  were  threefold. 
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1.  Since  passive  synthetic  aperture  can  be  viewed  as  a  scheme  that  converts 
temporal  gain  to  spatial  gain,  most  signals  of  interest  do  not  have 
sufficient  temporal  coherence  to  allow  a  long  spatially  coherent  aperture 
to  be  synthesized. 

2.  Since  past  algorithms  required  a  priori  knowledge  of  the  source 
frequency  in  order  to  compute  the  phase  correction  factor  [13],  the 
method  was  essentially  useless  in  any  bearing  estimation  problem  since 
Doppler  would  introduce  an  unknown  bias  on  the  frequency  observed  at 
the  receiver. 

3.  Since  synthetic  aperture  processing  essentially  converts  temporal  gain  to 
spatial  gain,  there  was  no  “new”  gain  to  be  achieved,  and  therefore,  no 
point  to  the  method. 

Recent  work  [12,13,60]  has  shown  that  there  can  be  realistic  conditions  in 
sonar  applications  under  which  all  of  these  objections  are  either  not  relevant 
or  do  not  constitute  serious  impediments  to  practical  applications  of  synthetic 
aperture  processing  in  operational  systems  [30],  Theoretical  discussions  have 
shown  [13]  that  the  above  three  arguments  are  valid  for  cases  that  include  the 
formation  of  synthetic  aperture  by  a  single  sensor  and  in  mediums  with 
isotropic  noise  characteristic.  However,  when  the  noise  characteristics  of  the 
received  signal  are  non-isotropic  and  the  receiving  array  includes  more  than 
one  sensor,  then  there  is  spatial  gain  available  from  passive  synthetic  aperture 
processing  and  this  has  been  discussed  analytically  in  [1 ,1 3,60],  Since  the 
medium  propagation  characteristics  of  sonar  and  ultrasound  system 
application  are  very  similar  and  non-isotropic,  synthetic  aperture  techniques 
that  are  successful  in  sonar  systems  can  be  successful  in  ultrasound 
applications  as  well.  Furthermore,  formation  of  a  synthetic  aperture  in  sonar 
systems  requires  the  movement  of  the  receiving  line  array.  In  ultrasound 
systems,  however,  the  formation  of  a  synthetic  aperture  can  be  derived  from 
the  sub-aperture  measurements  of  a  longer  receiving  array,  as  shown  in 
Figure  1,  which  eliminates  the  phase  errors  associated  with  the  movement  of 
the  receiving  array. 

A  summary  of  these  research  efforts  has  been  reported  in  a  special  issue  in 
the  IEEE  J.  Oceanic  Eng.  [13].  The  synthetic  aperture  processing  scheme 
that  has  been  used  in  broadband  sonar  applications  [1,60]  is  based  on  the 
Extended  Towed  Array  Measurements  (ETAM)  algorithm,  [12,30],  The  basic 
concept  of  this  algorithm  is  a  phase-correction  factor  that  is  used  to  combine 
coherently  successive  measurements  of  a  moving  receiving  line  array  to 
extend  the  effective  array  length. 

The  dark  shaded  sub-apertures  in  Figure  1  show  the  proposed  experimental 
implementation  of  the  ETAM  algorithm  for  ultrasound  applications  in  terms 
of  the  sub-aperture  line  array  size  and  sensor  positions  as  a  function  of  time 
and  space.  Between  two  successive  positions  of  the  32-sensor  sub-aperture 
there  are  a  number  of  sensor  pairs  of  space  samples  of  the  acoustic  field  that 
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have  the  same  spatial  information,  their  difference  being  a  phase  factor 
[1,12,13,60]  related  to  the  time  delay  these  measurements  were  taken.  The 
optimum  overlap  size,  which  is  related  to  the  variance  of  the  phase  correction 
estimates,  has  been  shown  [13]  to  be  equal  to  the  half  size  of  the  deployed 
sub-aperture.  For  the  particular  example  of  Figure  1,  the  spatial  overlap  size 
will  be  1 6-sensors.  Thus,  by  cross-correlating  the  16-sensor  pairs  of  the 
sensor  time  series  that  overlap,  the  desired  phase  correction  factor  is  derived, 
which  compensates  for  the  time  delay  between  these  measurements  and  the 
phase  fluctuations  caused  by  the  variability  and  non-isotropic  propagation 
characteristics  of  the  human  body;  this  is  called  the  overlap  correlator. 
Following  the  above,  the  key  parameters  in  the  ETAM  algorithm  is  the  time 
increment  rbetween  two  successive  sets  of  measurements.  This  may  be  the 
interval  of  0.3ms  between  two  active  ultrasound  transmissions.  Then,  the  total 
number  of  sets  of  measurements  required  by  the  32-scnsor  sub-aperture  to 
achieve  an  extended  aperture  size  equal  to  the  deployed  array  (i.e.  96-sensor 
array)  is  five. 

The  performance  characteristics  and  expectations  from  the  ETAM  algorithm 
have  been  evaluated  experimentally  for  sonar  applications  and  the  related 
results  have  been  reported  [12,13,30], 

Thus,  if  we  consider  the  ultrasound  system  in  Figure  1,  the  proposed  synthetic 
aperture  processing  will  coherently  synthesize  the  spatial  measurements 
derived  from  the  32-element  sub-apertures  of  the  ultrasound  receiving  array 
into  a  longer  aperture  equivalent  to  the  96-sensor  deployed  array  using  only  5 
sub-aperture  measurements  instead  of  64.  In  this  way,  the  required  hardware 
modifications  of  an  ultrasound  system  will  be  minimized  since  the  A/DC  will 
remain  the  same.  Moreover,  the  time  required  to  reconstruct  a  tomography 
image  will  be  reduced  from  the  current  20ms  time  interval  to  5x0.3  ms  =  1.5 
ms.  In  parallel  to  this  improvement,  there  will  be  an  increase  in  the  array  gain 
and  the  angular  resolution  of  the  ultrasound  beamforming  structure  by  5dBs. 


1.2.2  Adaptive  Beamforming  Structure  with  Near- 

Instantaneous  Convergence  for  Ultrasound  Systems 

The  adaptive  processing  schemes  proposed  in  this  investigation  are 
characterised  as  Dual-Use  technologies  and  have  been  tested  in  operational 
active  sonars  [60,61].  Real  data  results  have  shown  that  they  provide  array 
gain  improvements  for  signals  embedded  in  anisotropic  noise  fields,  similar 
to  those  in  the  human  body. 

Despite  the  geometric  differences  between  the  line  and  planar  arrays,  the 
underlying  beamforming  processes  for  these  arrays  are  time  delay 
beamforming  estimators,  which  are  basically  spatial  filters.  Flowever, 
optimum  beamforming  requires  the  beamforming  filter  coefficients  to  be 
chosen  based  on  the  covariance  matrix  of  the  received  data  by  the  N  -sensor 
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array  in  order  to  optimize  the  array  response  [60,61],  The  family  of 
algorithms  for  optimum  beamforming  that  use  the  characteristics  of  the  noise, 
are  called  Adaptive  Beamformers  [1,60].  In  this  investigation  we  consider  the 
implementation  in  ultrasound  systems  for  various  partially  adaptive  variants 
of  the  Linear  Constraint  Minimum  Variance  (LCMV)  method  and  a 
Generalized  Sidelobe  Canceller  (GSC)  adaptive  beamformer  [60,61]  that 
have  been  tested  successfully  in  real  time  active-passive  sonar  applications. 
Thus,  the  implementation  of  adaptive  schemes  in  ultrasound  systems  will  not 
be  restricted  into  one  method.  In  fact,  the  generic  concept  of  the  sub¬ 
aperture  multi-dimensional  array  to  be  considered  in  this  investigation  would 
allow  for  the  implementation  of  a  wide  variety  of  adaptive  schemes  in 
ultrasound  systems  [62],  As  for  the  implementation  of  adaptive  processing 
schemes  in  active  systems,  the  following  issues  need  to  be  addressed. 

For  ultrasound  applications  that  include  matched  filter  processing  for  Doppler 
shift  estimates,  the  outputs  of  the  adaptive  algorithms  are  required  to  provide 
coherent  beam  time  series  to  facilitate  the  post-processing.  This  means  that 
these  algorithms  should  exhibit  near-instantaneous  convergence  and  provide 
continuous  beam  time  series  that  have  sufficient  temporal  coherence  to 
correlate  with  the  reference  replica  in  matched  filter  processing  [1,60]. 


1.2.3  Advanced  Beamforming  Structure  for  Line  and  Planar 
Arrays  of  Ultrasound  Systems 

As  discussed  in  the  previous  section  (1.1.1),  deployment  of  planar  arrays  in 
ultrasound  medical  imaging  systems  has  been  gaining  increasing  popularity 
because  of  its  advantage  to  provide  real  3-D  images  of  organs  under  medical 
examination.  However,  operational  ultrasound  systems  deploying  planar 
arrays  are  not  yet  available.  Moreover,  if  we  consider  that  a  state-of-the-art 
phased  array  line  array  ultrasound  system  consists  of  128  sensors,  then  a 
planar  array  ultrasound  system  should  include  at  least  128x128  =  16,384 
sensors  in  order  to  achieve  the  angular  resolution  performance  of  a  line  array 
system  and  the  additional  3D  image  reconstruction  capability  provided  by  the 
elevation  beam  steering  of  a  planar  array.  Thus,  increased  angular  resolution 
in  azimuth  and  elevation  beam  steering  for  ultrasound  systems  means  larger 
sensor  arrays,  with  consequent  technical  and  higher  cost  implications.  As  it 
will  be  shown  in  this  report,  the  alternative  is  to  implement  synthetic  aperture 
and  adaptive  beam  processing  in  ultrasound  systems  that  deploy  a  planar 
array  with  1024  (64x16  or  32x32)  sensors,  which  consist  of  32  line  arrays 
with  32  sensors  each.  Then,  the  anticipated  array  gain  improvements  by  the 
adaptive  beamformer  would  be  equivalent  to  those  provided  by  a  96-sensor 
line  array  for  azimuth  beam-steering  and  a  96-sensor  vertical  line  array  for 
elevation  beam  steering  for  real  3D  ultrasound  imaging.  In  summary,  the 
array  gain  improvements  for  an  adaptive  1024-sensor  planar  array  will  be 
equivalent  to  those  that  could  be  provided  by  a  conventional  96x96  =  9216- 
sensor  planar  array.  This  is  because  for  line  arrays,  our  preliminary 
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quantitative  assessment  shows  that  the  image  resolution  improvements  of  the 
proposed  advanced  beamformers,  will  be  equivalent  to  a  three-time  longer 
physical  aperture. 


Figure  4.  Convential  3D  Beamforming  structure  for  2D  planar  array. 


To  address  this  requirement  and  to  allow  for  a  flexible  system  design  for  line 
and  or  planar  array  ultrasound  systems,  the  present  investigation  attempts  to 
develop  an  advanced  beamforming  structure  applicable  to  both  line  and 
planar  arrays.  The  proposed  approach  is  conceptually  similar  to  previous 
investigations  [7,26,61]  carried  out  for  sonar  system  applications  by  the 
authors  of  this  report. 

Briefly,  these  previous  studies  have  shown  that  a  planar  array  beamformer 
can  be  decomposed  into  tw'o  line  array  beamforming  steps.  For  example,  let 
us  consider  the  discrete  planar  array  in  Figure  4,  w'ith  N  sensors  where 
N  =  NM  and  M,  N  are  the  number  of  sensors  along  x-axis  and  y-axis, 
respectively.  The  first  step  includes  an  M-sensor  line-array  beamforming 
along  the  x-axis  of  the  planar  array  providing  beam  time  series  with  azimuth 
beam  steering.  The  line  array  beamformer  will  be  repeated  N-  times  to  get  the 
vector  that  includes  the  beam  times  series  By(f,6 5)  ,  w'here  the  index 

y  =  —  1  is  along  the  axis  y.  For  a  given  azimuth  beam  steering,  the 

second  step  includes  line  array  beamforming  along  y-axis  by  treating  the 
vector  Br{f,Os)  as  the  input  signal  for  the  line  array  beamformer  to  get  the 

output  B(  f,  9 s ,  <ps ) ,  which  is  the  output  of  a  planar  array  beamformer  for  a 
given  azimuth  9S  and  elevation  beamsteering  <j)s ,  respectively. 

Figure  4  show's  the  involved  steps  of  decomposing  the  2-D  planar  array 
beamformer  into  two  steps  of  line-array  beamformers.  The  decomposition  of 
the  planar  array  beamformer  into  these  tw'o  line-array  beamforming  steps 
leads  to  an  efficient  implementation  based  on  the  following  tw'o  factors. 
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First,  the  number  of  the  involved  sensors  for  each  of  these  line  array 
beamformers  is  much  less  than  the  total  number  of  sensors,  X  of  the  planar 
array.  This  type  of  decomposition  process  for  the  2-D  beamformer  eliminates 
the  need  for  very  large  memory  and  CPU's  with  very  high  throughput 
requirements  in  one  board  for  real  time  system  applications.  Secondly,  all 
these  line  array  beamformers  can  be  executed  in  parallel,  which  allows  their 
implementation  in  much  simpler  parallel  architectures  with  simpler  CPU's, 
which  is  a  practical  requirement  for  real  time  system  application.  Besides  the 
advantage  of  the  efficient  implementation,  the  proposed  decomposition 
approach  makes  the  application  of  the  spatial  window  much  simpler  to  be 
incorporated.  Furthermore,  both  the  synthetic  aperture  and  adaptive 
processing  schemes  of  Sections  2.3.1  and  2.3.2  that  will  be  developed  for  line 
arrays,  will  be  directly  applicable  to  planar  arrays  through  the  proposed 
decomposition  of  the  planar  array  beamforming  process  into  two  line  array 
beamforming  steps. 
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2.  Theoretical  Remarks 


Active  ultrasound  operations  involve  the  transmission  of  well  defined  ultrasound  signals, 
called  replicas,  which  illuminate  the  human  body  medium.  The  reflected  ultrasound  energy 
from  an  extended  target  or  body  organ  provides  the  array  receiver  with  a  basis  for  detection 
and  estimation.  The  image  reconstruction  process  is  a  passive  ultrasound  bcamforming 
operation. 

All  the  above  active-passive  beamforming  issues  have  been  discussed  in  several  review 
articles  [1-6]  that  form  a  good  basis  for  interested  readers  to  become  familiar  with  “main 
stream”  beamforming  signal  processing  developments.  Therefore,  discussions  of  issues  of 
conventional  signal  processing,  detection,  estimation  and  influence  of  medium  on  ultrasound 
system  performance  are  beyond  the  scope  of  this  report.  Only  a  very  brief  overview  of  the 
above  issues  will  be  highlighted  in  this  section  in  order  to  define  the  basic  terminology 
required  for  the  presentation  of  the  main  theme  of  the  present  article.  Let  us  start  with  a 
basic  system  model  that  reflects  the  interrelationships  between  the  extended  target  or  the 
human  body  (medium)  and  the  receiving  sensor  array  of  an  ultrasound  system. 

A  schematic  diagram  of  this  basic  system  is  shown  in  Figure  5,  where  array  signal 
processing  is  shown  to  be  two-dimensional  [1,5,10,12,18]  in  the  sense  that  it  involves  both 
temporal  and  spatial  spectral  analysis.  The  temporal  processing  provides  spectral 
characteristics  that  are  used  for  classification  and  the  spatial  processing  provides  estimates 
of  the  directional  characteristics,  (i.e.  bearing  and  depth),  of  a  detected  signal.  Thus,  Space- 
Time  Processing  is  the  fundamental  processing  concept  in  ultrasound  systems  and  it  will  be 
the  subject  of  our  discussion  in  the  next  section. 

2.1  Space-Time  Processing 

For  geometrical  simplicity  and  without  any  loss  of  generality,  we  consider  here  a 
combination  of  N  equally  spaced  sensors  in  a  linear  array,  which  may  form  a  phased  array 
system  that  can  be  used  to  estimate  the  directional  properties  of  echoes  and  acoustic  signals. 
As  shown  in  Fig.  5,  a  direct  analogy  between  sampling  in  space  and  sampling  in  time  is  a 
natural  extension  of  the  sampling  theory  in  space-time  signal  representation  and  this  type  of 
space-time  sampling  is  the  basis  in  array  design  that  provides  a  description  of  an  array 
system  response.  When  the  sensors  are  arbitrarily  distributed,  each  element  will  have  an 
added  degree  of  freedom,  which  is  its  position  along  the  axis  of  the  array. 

This  is  analogous  to  non-uniform  temporal  sampling  of  a  signal.  In  this  report  we  extend  our 
discussion  to  multi-dimensional  array  systems. 

Sources  of  sound  that  are  of  interest  in  ultrasound  system  applications  provide  harmonic 
narrowband  and  broadband  signals  that  satisfy  the  wave  equation  [2,10],  Furthermore,  their 
solutions  have  the  property  that  their  associated  temporal-spatial  characteristics  are 
separable  [10],  Therefore,  measurements  of  the  pressure  field  z(?,t)  which  is  excited  by 
acoustic  source  signals,  provide  the  spatial-temporal  output  response,  designated  by  x(r,t ) 
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of  the  measurement  system.  The  vector  r  refers  to  the  source-sensor  relative  position  and  t 
is  the  time. 


Figure  5.  A  model  of  space-time  signal  processing.  It  shows  that  ultrasound  signal  processing  is  two 
dimensional  in  the  sense  that  it  involves  both  temporal  and  spatial  spectral  analysis.  The  temporal 
processing  provides  characteristics  for  target  classification  and  the  spatial  processing  provides 
estimates  of  the  directional  characteristics  (bearing,  range-depth)  of  detected  echoes  (active  case)  or 

signals  of  interest  (passive  case). 


The  output  response  is  the  convolution  of  z(r>0  with  the  line  array  system 

response  h{r,t)  [ 1 0,30] 


x(r,t )  =  z(r,t)®h(r,t) 


(4) 


where  0  refers  to  convolution.  Since  z(r»0  is  defined  at  the  input  of  the  receiver,  it  is  the 
convolution  of  the  source's  characteristics  0  with  the  underwater  medium's  response 


z(r,t)  =  y(r,f)®  yAjj) 

Fourier  transformation  of  Equation  (3)  provides: 

X(Q),k)=  ^{a),k)^{(o,k)\  H(co,k) 


(5) 

(6) 
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where,  Q),k  are  the  frequency  and  wavenumber  parameters  of  the  temporal  and  spatial 
spectrums  of  the  transform  functions  in  Equations  (4)  &  (5).  Signal  processing,  in  terms  of 
beamforming  operations,  of  the  receiver's  output  x(r,t )  ,  provides  estimates  of  the  source 
bearing  and  possibly  of  the  source  range.  This  is  a  well-understood  concept  of  the  forward 
problem,  which  is  concerned  with  determining  the  parameters  of  the  received  signal  x(r,t) 
given  that  we  have  information  about  the  other  two  functions  z(r,t )  and  [5],  The 

inverse  problem  is  concerned  with  determining  the  parameters  of  the  impulse  response  of 
the  medium  y/(r,t)  by  extracting  information  from  the  received  signal  x(r,t)  assuming 
that  the  function  h(r,t )  is  known  [5],  The  ultrasound  problems,  however,  are  quite 
complex  and  include  both  forward  and  inverse  problem  operations.  In  particular,  detection, 
estimation  and  tracking-localization  processes  of  ultrasound  systems  are  typical  examples  of 
the  forward  problem,  while  target  classification  for  passive-active  diagnostic  ultrasound 
imaging  are  typical  examples  of  the  inverse  problem.  In  general,  the  inverse  problem  is  a 
computationally  very  costly  operation  and  typical  examples  in  acoustic  signal  processing  are 
seismic  deconvolution  and  acoustic  tomography. 


2.2  Definition  of  Basic  Parameters 

This  section  outlines  the  context  in  which  the  ultrasound  problem  can  be  viewed  in  terms  of 
simple  models  of  acoustic  signals  and  noise  fields.  The  signal  processing  concepts  that  are 
discussed  in  this  report  have  been  included  in  sonar  and  radar  investigations  with  sensor 
arrays  having  circular,  planar,  cylindrical  and  spherical  geometric  configurations 
[7,25,26,28].  Thus,  we  consider  a  multi-dimensional  array  of  equally  spaced  sensors  with 
spacing  8.  The  output  of  the  «th  sensor  is  a  time  series  denoted  by  xn(tj)  ,  where  Ms 

)  are  the  time  samples  for  each  sensor  time  series.  *  denotes  complex  conjugate  transposition 
—  * 

so  that  X  is  the  row  vector  of  the  received  N  -  sensor  time  series  {xn(t  -J  ,n=],2,...,  X  }. 


Then  xn(tj)  =  sn(tj)  +  £n(tj)  ,  where  sn(tj)  ,  £n(tj)  are  the  signal  and  noise  components  in  the 
received  sensor  time  series.  5  ,  £  denote  the  column  vectors  of  the  signal  and  noise 
components  of  the  vector  x  of  the  sensor  outputs  (i.e.  x  =  .v+f  )• 


Ms 

Xn(f)  =  T  xn  (/, )  exp(- j2jfti ) 

is  the  Fourier  transform  of  xn(tj )  at  the  signal  with 
frequency/  c  —fX  is  the  speed  of  sound  in  the  underwater,  or  human-body  medium  and  X 
is  the  wavelength  of  the  frequency/  5  =  £’{5'  5*}  is  the  spatial  correlation  matrix  of 
the  signal  vector  $  »  whose  nth  element  is  expressed  by, 


(*,■)  =  *„[ >,-+*■„  (0,0)]  , 


(7) 


E{...}  denotes  expectation  and  T n  (0 ,0)  is  the  time  delay  between  the  (n-l)st  and  the  nth 
sensor  of  the  array  for  an  incoming  plane  wave  with  direction  of  propagation  of  azimuth 
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angle  9  and  an  elevation  angle  <p,  as  depicted  in  Figure  3.  In  the  frequency  domain,  the 
spatial  correlation  matrix  S  for  the  plane  wave  signal  sn(tj)  is  defined  by: 


S(Jt  ,0,0)  =  As  <Jt  )D  (/, ,  0,  0)D  *  <Jt  ,6,0) , 


(8) 


where  As{ft)  is  the  power  spectral  density  of  s(tj)  for  the  ith  frequency  bin;  and 
D{f,6,(p)  is  the  steering  vector  with  the  n*  term  being  denoted  by  dn(f ,9,0) .  Then 
matrix  S(fj9,(f))  has  as  its  nth  row  and  with  column  defined  by,  Snm(f\,  9,(f))  - 
A s (fi)dn (fi 9 <j>) d *m (fi d,  <ft) .  Moreover,  R(f, )  is  the  spatial  correlation  matrix  of  received 

sensor  time  series  with  elements,  R„m(  f, dnm )  .  Re  (fi  )  =  <7 n  (/]  )R£  (fl )  is  the  spatial 
correlation  matrix  of  the  noise  for  the  ith  frequency  bin  with  o\  (/,.)  being  the  power 
spectral  density  of  the  noise,  £„(tj).  In  what  is  considered  as  an  estimation  procedure  in  this 
report,  the  associated  problem  of  detection  is  defined  in  the  classical  sense  as  a  hypothesis 
test  that  provides  a  detection  probability  and  a  probability  of  false  alarm  [31-33].  This 
choice  of  definition  is  based  on  the  standard  CFAR  (constant  false  alarm  rate)  processor, 
which  is  based  on  the  Neyman-Pearson  criterion  [31].  The  CFAR  processor  provides  an 
estimate  of  the  ambient  noise  or  clutter  level  so  that  the  threshold  can  be  varied  dynamically 
to  stabilize  the  false  alarm  rate.  Ambient  noise  estimates  for  the  CFAR  processor  are 
provided  mainly  by  noise  normalization  techniques  [34]  that  account  for  the  slowly  varying 
changes  in  the  background  noise  or  clutter.  The  above  estimates  of  the  ambient  noise  are 
based  upon  the  average  value  of  the  received  signal,  the  desired  probability  of  detection  and 
probability  of  false  alarms. 

At  this  point,  a  brief  discussion  on  the  fundamentals  of  detection  and  estimation  process  is 
required  in  order  to  address  implementation  issues  of  signal  processing  schemes  in 
ultrasound  systems. 


2.3  Detection  and  Estimation 

In  passive  systems,  in  general,  we  do  not  have  the  a  priori  probabilities  associated  with  the 
hypothesis  H]  that  the  signal  is  assumed  present  and  the  null  hypothes  is  Ho  that  the 
received  time  series  consists  only  of  noise.  As  a  result,  costs  can  not  be  assigned  to  the 
possible  outcomes  of  the  experiment.  In  this  case,  the  Neyman-Pearson  (N-P)  Criterion  [31] 
is  applied  because  it  requires  only  a  knowledge  of  the  signal's  and  noise's  probability  density 
functions  (pdf). 

Let  xn=i(tj)  ,  ( i=l,...,M )  denote  the  received  vector  signal  by  a  single  sensor.  Then  for 
hypothesis  Hi,  which  assumes  that  the  signal  is  present,  we  have: 

Hi:  xn=lOi)  ~  sn=l9i)  7£n=](tj)  , 
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where  sn-j(tj)  and  £n=i(tj)  are  the  signal  and  noise  vector  components  in  the  received  signal 
and  pj(x)  is  the  pdf  of  the  received  signal  xn=i(tj)  given  that  Hj  is  true.  Similarly,  for 
hypothesis  Hq  : 


H():  xn=lOi)  £n=l  Oi) 


and  po(x)  is  the  pdf  of  the  received  signal  given  that  Hq  is  true.  The  N-P  criterion  requires 
maximization  of  probability  of  detection  for  a  given  probability  of  false  alarm.  So,  there 
exists  a  non-negative  number  77  such  that  if  hypothesis  Hj  is  chosen  then 


A(x)  = 


,w 


>n 


(9) 


which  is  the  likelihood  ratio.  By  using  the  analytic  expressions  for  po(x)  (the  pdf  for  Hq) 
and  pi(x)  (the  pdf  for/// )  in  Eq.  (9)  and  by  taking  the  In  [X(x)]  ,  we  have  [31], 


Ar  -  In  [A (x )]  =  s*R£x 


(10) 


where,  Xt  is  the  log  likelihood  ratio  and  R  £  is  the  covariance  matrix  of  the  noise  vector, 

as  defined  in  the  previous  section  2.2.  For  the  case  of  white  noise  with  R£=(*>J  and/ 
the  unit  matrix,  the  test  statistic  in  expression  (10)  is  simplified  into  a  simple  correlation 
receiver  (or  replica  correlator) 


Ar  =  s'  ®  X  .  (11) 

For  the  case  of  anisotropic  noise,  however,  an  optimum  detector  should  include  the 
correlation  properties  of  the  noise  in  the  correlation  receiver  as  this  is  defined  in  Eq.  (10). 

For  plane  wave  arrivals  that  are  observed  by  an  jV -sensor  array  receiver  the  test  statistics 
are  [31]: 


K  =  £  (/,) •*!(/,)■  [s(y , <m)  +  r!c (f )]-’  x(f:) 

(=i  ,v 1 H 

where,  the  above  statistics  are  for  the  frequency  domain  with  parameters  defined  in  Eqs.  (7), 
(8)  in  the  previous  section  2.2.  Then,  for  the  case  of  an  array  of  sensors  receiving  plane 
wave  signals,  the  log  likelihood  ratio  Trin  Eq.  (12)  is  expressed  by  the  following  equation, 


18 


DR  DC  Toronto  TR  2002-058 


which  is  the  result  of  simple  matrix  manipulations  based  on  the  frequency  domain 
expressions  (7),  (8)  and  their  parameter  definitions  presented  in  Section  2.2.  Thus, 


M 


-1 


k  =  y,  rw,) 


i=i 


(13) 


where  [31], 

v2(f,)  = 


4  (/,)/<*;(/,) 


1  +  A,  if,  )D‘  (/ ,  </>,  e  y;'  (/ )!)(/)/  a]  (/, ) 


Eq.  (13)  can  be  written  also  as  follows, 


(14) 


J.-I 

/= 1 


N 


(15) 


This  last  expression  (15)  of  the  log  likelihood  ratio  indicates  that  an  optimum  detector  in  this 
case  requires  the  filtering  of  each  one  of  the  N  -sensor  received  time  series  X„  (fi)  with  a  set 
of  filters  being  the  elements  of  the  vector, 


(16) 


Then,  the  summation  of  the  filtered  sensor  outputs  in  frequency  domain  according  to  Eq. 
(16)  provides  the  test  statistics  for  optimum  detection.  For  the  simple  case  of  white  noise 

=  &  n  I  and  for  a  line  array  receiver,  the  filtering  operation  in  (16)  indicates  plane 
wave  conventional  beamforming  in  the  frequency  domain, 


N 


(17) 


where,  —  £  /(l  +  Ng ) 

Q  =  Altai. 


is  a  scalar,  which  is  a  function  of  the  signal-to-noise  ratio, 


For  the  case  of  narrowband  signals  embedded  in  spatially  and/or  temporally  correlated 
noise,  it  has  been  shown  [13]  that  the  deployment  of  very  long  arrays  or  application  of 
acoustic  synthetic  aperture  will  provide  sufficient  array  gain  and  will  achieve  optimum 
detection  and  estimation  for  the  parameters  of  interest. 
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For  the  general  case  of  broadband  and  narrowband  signals  embedded  in  a  spatially 
anisotropic  and  temporally  correlated  noise  field,  expression  (17)  indicates  that  the  filtering 
operation  for  optimum  detection  and  estimation  requires  adaptation  of  the  ultrasound  signal 
processing  according  to  the  human  body’s  noise  characteristics,  respectively.  The  family  of 
algorithms  for  optimum  beamforming  that  use  the  characteristics  of  the  noise,  are  called 
Adaptive  Beamformers  [3,17-20,22,23];  and  a  detailed  definition  of  an  adaptation  process 

requires  knowledge  of  the  correlated  noise's  covariance  matrix  (/,  )  .  However,  if  the 
required  knowledge  of  the  noise’s  characteristics  is  inaccurate,  the  performance  of  the 
optimum  beamformer  will  degrade  dramatically  [1 8,23],  As  an  example,  the  case  of 
cancellation  of  the  desired  signal  is  often  typical  and  significant  in  adaptive  beamforming 
applications  [1 8,24].  This  suggests  that  the  implementation  of  useful  adaptive  beamformers 
in  real  time  operational  systems  is  not  a  trivial  task.  The  existence  of  numerous  articles  on 
adaptive  beamforming  suggests  the  dimensions  of  the  difficulties  associated  wdth  this  type  of 
implementation.  In  order  to  minimize  the  generic  nature  of  the  problems  associated  with 
adaptive  beamforming  the  concept  of  partially  adaptive  beamformer  design  w'as  introduced. 
This  concept  reduces  the  degrees  of  freedom,  which  results  in  lowering  the  computational 
requirements  and  often  improving  the  adaptive  response  time  [1 7,1 8],  However,  the  penalty 
associated  with  the  reduction  of  the  degrees  of  freedom  in  partially  adaptive  beamformers  is 
that  they  cannot  converge  to  the  same  optimum  solution  as  the  fully  adaptive  beamformer. 

Although  a  review  of  the  various  adaptive  beamformers  would  seem  relevant  at  this  point, 
we  believe  that  this  is  not  necessary  since  there  are  excellent  review  articles  [3,17,18,21]  that 
summarize  the  points  that  have  been  considered  for  the  material  of  this  report.  There  are 
two  main  families  of  adaptive  beamformers,  the  Generalized  Sidelobe  Cancellers  (GSC) 
[44,45]  and  the  Linearly  Constrained  Minimum  Variance  Beamformers  (LCMV)  [18].  A 
special  case  of  the  LCMV  is  Capon's  Maximum  Likelihood  Method  [22],  which  is  called 
Minimum  Variance  Distortionless  Response  (MVDR)  [17,18,22,23,38,39],  This  algorithm 
has  proven  to  be  one  of  the  more  robust  of  the  adaptive  array  beamformers  and  it  has  been 
used  by  numerous  researchers  as  a  basis  to  derive  other  variants  of  MVDR  [18],  In  this 
report  we  will  address  implementation  issues  for  various  partially  adaptive  variants  of  the 
MVDR  and  a  GSC  adaptive  beamformer  [1],  w'hich  are  discussed  in  Section  4.2. 

In  summary,  the  classical  estimation  problem  assumes  that  the  a  priori  probability  of  the 
signal's  presence  p(Hj)  is  unity  [31-33],  However,  if  the  signal's  parameters  are  not 
known  a  priori  and  p  (H/ )  is  knowm  to  be  less  than  unity,  then  a  series  of  detection 
decisions  over  an  exhaustive  set  of  source  parameters  constitutes  a  detection  procedure, 
where  the  results  incidentally  provide  an  estimation  of  source's  parameters.  As  an  example, 
we  consider  the  case  of  a  matched  fdter,  which  is  used  in  a  sequential  manner  by  applying  a 
series  of  matched  filter  detection  statistics  to  estimate  the  range  and  speed  of  the  target, 
which  are  not  known  a  priori .  This  type  of  estimation  procedure  is  not  optimal  since  it  does 
not  constitute  an  appropriate  form  of  Bayesian  minimum  variance  or  minimum  mean  square 
error  procedure. 

Thus,  the  problem  of  detection  [31-33]  is  much  simpler  than  the  problem  of  estimating  one 
or  more  parameters  of  a  detected  signal.  Classical  decision  theory  [31-33,]  treats  signal 
detection  and  signal  estimation  as  separate  and  distinct  operations.  A  detection  decision  as 
to  the  presence  or  absence  of  the  signal  is  regarded  as  taking  place  independently  of  any 
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signal  parameter  or  waveform  estimation  that  may  be  indicated  as  the  result  of  detection 
decision.  However,  interest  in  joint  or  simultaneous  detection  and  estimation  of  signals 
arises  frequently.  Middleton  and  Esposito  [46]  have  formulated  the  problem  of  simultaneous 
optimum  detection  and  estimation  of  signals  in  noise  by  viewing  estimation  as  a  generalized 
detection  process.  Practical  considerations,  however,  require  different  cost  functions  for 
each  process  [46],  As  a  result,  it  is  more  effective  to  retain  the  usual  distinction  between 
detection  and  estimation. 

Estimation,  in  ultrasound  systems,  includes  both  the  temporal  and  spatial  structure  of  an 
observed  signal  field.  For  active  systems,  correlation  processing  and  Doppler  (for  moving 
target  indications)  are  major  concerns  that  define  the  critical  distinction  between  these  two 
approaches  (i.e.  passive,  active )  to  ultrasound  processing.  In  this  report,  we  restrict  our 
discussion  only  to  topics  related  to  spatial  signal  processing  for  estimating  signal 
parameters,  which  is  associated  with  the  image  reconstruction  process  of  ultrasound 
systems.  However,  spatial  signal  processing  has  a  direct  representation  that  is  analogous  to 
the  frequency-domain  representation  of  temporal  signals.  Therefore,  the  spatial  signal 
processing  concepts  discussed  here  have  direct  applications  to  temporal  spectral  analysis. 


2.4  Cramer-Rao  Lower  Bound  (CRLB)  Analysis 

Typically,  the  performance  of  an  estimator  is  represented  as  the  variance  in  the  estimated 
parameters.  Theoretical  bounds  associated  with  this  performance  analysis  are  specified  by 
the  Cramer-Rao  bound  [31-33]  and  that  has  led  to  major  research  efforts  by  the  sonar  signal 
processing  communi  ty  in  order  to  define  the  idea  of  an  optimum  processor  for  discrete 
sensor  arrays  [12,16,  56-60].  If  the  a  priori  probability  of  detection  is  close  to  unity  then  the 
minimum  variance  achievable  by  any  unbiased  estimator  is  provided  by  the  Cramer-Rao 
Lower  Bound  (CRLB)  [3 1,32,46], 

More  specifically,  let  us  consider  that  the  received  signal  by  the  nth  sensor  of  a  receiving 
array  is  expressed  by, 


xn0i)  ~  sn(ij)  F  ^nOi) 


(18) 


where,  s„  (/,  >  © )  —  s „  I X  +  (0, 0)]  s  defines  the  received  signal  model  with 

T  n  (0  ,(p  )  being  the  time  delay  between  the  (n-l)st  and  the  nth  sensor  of  the  array  for  an 
incoming  plane  wave  with  direction  of  propagation  of  azimuth  angle  and  an  elevation  angle 

0 ,  as  depicted  in  Figure  3.  The  vector  0  ,  includes  all  the  unknown  parameters  considered 

2 

in  relation  (18).  Let  O  e_  denote  the  variance  of  an  unbiased  estimate  of  an  unknown 
parameter  in  the  vector  0  .  The  Cramer-Rao  [31-33]  bound  states  that  the  best  unbiased 
estimate  0  of  the  parameter  vector  0  has  the  covariance  matrix 


DRDC  Toronto  TR  2002-058 


21 


cov©  >  y(©) ' 


(19) 


where  J  is  the  Fisher  information  matrix  whose  elements  are: 


J,=~E\ 


'  if  In  p(x 

e)" 

mimi 

V 

J 

(20) 


In  Eq.  (20), 


P  X 


0 


is  the  probability  density  function  (pdf)  governing  the  observations: 


X  —  [Xj  (/,.  ),x2  (/,  ),x3  (/,  ),...xA;  (/,.  )J  5  for  each  of  the  N  and  Ms  independent  spatial  and 
temporal  samples  respectively  that  are  described  by  the  model  in  Eq.  (18).  The  variance  of 
the  unbiased  estimates  0  has  a  lower  bound  (called  the  CRLB),  which  is  given  by  the 
diagonal  elements  of  expression  (19).  This  CRLB  is  used  as  standard  of  performance  and 
provides  a  good  measure  for  the  performance  of  a  signal-processing  algorithms  which  gives 

unbiased  estimates  0  for  the  parameter  vector  0  .  In  this  case,  if  there  exists  a  signal 
processor  to  achieve  the  CRLB,  it  will  be  the  maximum-likelihood  estimation  (MLE) 
technique.  The  above  requirement  associated  with  the  a  priori  probability  of  detection  is 
very  essential  because  if  it  is  less  than  one,  then  the  estimation  is  biased  and  the  theoretical 
CRLBs  do  not  apply.  This  general  framework  of  optimality  is  very  essential  in  order  to 
account  for  Middleton's  [32]  warning  that  a  system  optimized  for  the  one  function  (detection 
or  estimation)  may  not  be  necessarily  optimized  for  the  other.  For  a  given  model  describing 
the  received  signal  by  a  sonar  or  ultrasound  system,  the  CRLB  analysis  can  be  used  as  a  tool 
to  define  the  information  inherent  in  a  sonar  system.  This  is  an  important  step  related  to  the 
development  of  the  signal  processing  concept  for  a  sonar  system  as  well  as  in  defining  the 
optimum  sensor  configuration  arrangement  under  which  we  can  achieve,  in  terms  of  system 
performance,  the  optimum  estimation  of  signal  parameters  of  our  interest.  This  approach  has 
been  applied  successfully  to  various  studies  related  to  the  present  development  [12,15,56- 
60]. 


As  an  example,  let  us  consider  the  simplest  problem  of  one  source  with  bearing  0,  being 

the  unknown  parameter.  Following  relation  (20),  the  results  of  the  variance  <7^  jn  the 
bearing  estimates  are, 


°e.  = 


2  y/N 


B. 


V 


n  sin  6 


i  J 


(21) 


where,  Iff  —  M  SA^  /  O' N  ,  the  parameter  Bw  —Xl(N  —  X)8  gives  the 
beamwidth  of  the  physical  aperture  that  defines  the  angular  resolution  associated  with  the 
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estimates  of  6X .  The  signal  to  noise  ratio  (SNR)  at  the  sensor  level  is 
SNR  =  10xlog10(yO  or 

SNR  =  20  x  log  10  (Ax  / <7, )  +  10  x  log  J0(MS)  (22) 

2 

It  is  obvious  from  the  above  relations  (21)  and  (22)  that  the  variance  of  the  bearing  C7  e  can 
get  smaller  when  the  observation  period  T  =  MS/  fs  becomes  long  and  the  receiving 
array  size,  L  =  ( N  —  l)/l  gets  very  long  . 

The  next  question  needed  to  be  addressed  is  about  the  unbiased  estimator  that  can  exploit 
this  available  information  and  provide  results  asymptotically  reaching  the  CRLBs.  For  each 
estimator  it  is  well  known  that  there  is  a  range  of  Signal-to-Noise  Ratio  (SNR)  in  which  the 
variance  of  the  estimates  rises  very  rapidly  as  SNR  decreases.  This  effect,  which  is  called  the 
threshold  effect  of  the  estimator,  determines  the  range  of  SNR  of  the  received  signals  for 
which  the  parameter  estimates  can  be  accepted.  In  passive  sonar  systems  the  SNR  of  signals 
of  interest  are  often  quite  low  and  probably  below  the  threshold  value  of  an  estimator.  In  this 
case,  high  frequency  resolution  in  both  time  and  spatial  domains  for  the  parameter 
estimation  of  narrowband  signals  is  required.  In  other  words,  the  threshold  effect  of  an 
estimator  determines  the  frequency  resolution  for  processing  and  the  size  of  the  array 
receivers  required  in  order  to  detect  and  estimate  signals  of  interest  that  have  very  low  SNR 
[12,14,53,61,62].  The  CRLB  analysis  has  been  used  in  many  studies  to  evaluate  and 
compare  the  performance  of  the  various  non-conventional  processing  schemes  [17,18,55] 
that  have  been  considered  for  implementation  in  the  generic  beamforming  structure  to  be 
discussed  in  Section  4.1.  In  general,  array  signal  processing  includes  a  large  number  of 
algorithms  for  a  variety  of  systems  that  are  quite  diverse  in  concept.  There  is  a  basic  point 
that  is  common  in  all  of  them,  however,  and  this  is  the  beamforming  process,  which  we  are 
going  to  examine  in  Section  3. 
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3.  Optimum  Estimators  for  Array  Signal  Processing 


An  one  dimensional  (1-D)  device  such  as  a  line  sensor  array  satisfies  the  basic  requirements 
of  a  spatial  filter.  It  provides  direction  discrimination,  at  least  in  a  limited  sense,  and  a 
signal-to-noise  ratio  improvement  relative  to  an  omni-directional  sensor.  Because  of  the 
simplified  mathematics  and  reduced  number  of  the  involved  sensors,  relative  to  multi¬ 
dimensional  arrays,  most  of  the  researchers  have  focused  on  the  investigation  of  the  line 
sensor  arrays  in  system  applications  [1-6].  Furthermore,  implementation  issues  of  synthetic 
aperture  and  adaptive  techniques  in  real  time  systems  have  been  extensively  investigated  for 
line  arrays  as  well  [1,5,6,12,17,19,20],  However,  the  configuration  of  the  array  depends  on 
the  purpose  for  which  it  is  to  be  designed.  For  example,  if  a  wide  range  of  horizontal  angles 
is  to  be  observed,  a  circular  configuration  may  be  used,  that  provides  beam  characteristics 
that  are  independent  of  the  direction  of  steering.  Vertical  direction  may  be  added  by  moving 
into  cylindrical  configuration  [8],  In  a  more  general  case,  where  both  vertical  and 
horizontal  steering  is  to  be  required  and  where  a  large  range  of  angles  is  to  be  covered,  a 
spherically  symmetric  array  would  be  desirable  [9],  In  modern  ultrasound  imaging  systems, 
planar  arrays  are  required  to  reconstruct  real-time  3-D  images.  However,  the  huge 
computational  load  required  for  multi-dimensional  conventional  and  adaptive  beamformers 
makes  the  applications  of  these  2-D  &  3-D  arrays  in  real-time  systems  non  feasible. 

Furthermore,  for  modern  ultrasound  systems,  it  has  become  a  necessity  these  days  that  all 
possible  active  and  passive  modes  of  operation  should  be  exploited  under  an  integrated 
processing  structure  that  reduces  redundancy  and  provides  cost  effective  real  time  system 
solutions  [6].  Similarly,  the  implementation  of  computationally  intensive  data  adaptive 
techniques  in  real  time  systems  is  also  an  issue  of  equal  practical  importance.  Flowever, 
when  these  systems  include  multi-dimensional  (2-D,  3-D)  arrays  with  hundreds  of  sensors, 
then  the  associated  beamforming  process  requires  very  large  memory  and  very  intensive 
throughput  characteristics,  something  that  makes  its  implementation  in  real  time  systems  a 
very  expensive  and  difficult  task. 

To  counter  this  implementation  problem,  the  present  report  introduces  a  generic  approach  of 
implementing  conventional  beamforming  processing  schemes  with  integrated  passive  and 
active  modes  of  operations  in  systems  that  may  include,  planar,  cylindrical  or  spherical 
arrays  [25-28],  This  approach  decomposes  the  2-D  and  3-D  beamforming  process  into  sets 
of  line  and/or  circular  array  beamformers.  Because  of  the  decomposition  process,  the  fully 
multi-dimensional  beamformer  can  now  be  divided  into  sub-sets  of  coherent  processes  that 
can  be  implemented  in  small  size  CPU’s  that  can  be  integrated  under  the  parallel 
configuration  of  existing  computing  architectures.  Furthermore,  application  of  spatial 
shading  for  multidimensional  beamformers  to  control  side-lobe  structures  can  now  be  easily 
incorporated.  This  is  because  the  problem  of  spatial  shading  for  line  arrays  has  been 
investigated  thoroughly  [36]  and  the  associated  results  can  be  integrated  into  a  circular  and  a 
multi-dimensional  beamformer,  which  can  be  decomposed  now  into  coherent  sub-sets  of 
line  and  or  circular  beamformers  of  the  proposed  generic  processing  structure. 

As  a  result  of  the  decomposition  process,  provided  by  the  generic  processing  structure,  the 
implementation  effort  for  adaptive  schemes  is  reduced  to  implementing  adaptive  processes 
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in  line  and  circular  arrays.  Thus,  a  multi-dimensional  adaptive  beamformer  can  now  be 
divided  into  two  coherent  modular  steps  which  lead  to  efficient  system  oriented 
implementations.  In  summary,  the  proposed  approach  demonstrates  that  the  incorporation  of 
adaptive  schemes  with  near-instantaneous  convergence  in  multi-dimensional  arrays  is 
feasible  [7,25-28], 

At  this  point  it  is  important  to  note  that  the  proposed  decomposition  process  of  2-D  and  3-D 
conventional  beamformers  into  sets  of  line  and/or  circular  array  beamformers  is  an  old 
concept  that  has  been  exploited  over  the  years  by  sonar  system  designers.  Thus,  references 
on  this  subject  may  exist  in  reports  of  Navy  laboratories  and  industrial  institutes  that  are  not 
always  readily  available  and  the  authors’  of  this  document  are  not  aware  of  any  type  of 
reports  in  this  area.  Previous  efforts  attempted  to  address  practical  implementation  issues 
and  had  been  focused  on  cylindrical  arrays.  As  an  example,  a  cylindrical  array  beamformer 
is  decomposed  into  time-delay  line  array  beamformers  providing  beams  along  elevation 
angles  of  the  cylindrical  array.  These  are  called  staves.  Then,  the  beam  time  series  associated 
with  a  particular  elevation  steering  of  interest  are  provided  at  the  input  of  a  circular  array 
beamformer. 

In  this  report  the  attempt  is  to  provide  a  higher  degree  of  development  than  the  one 
discussed  above  for  cylindrical  arrays.  The  task  is  to  develop  a  generic  processing  structure 
that  integrates  the  decomposition  process  of  multi-dimensional  planar,  cylindrical  and 
spherical  array  beamformers  into  line  and  or  circular  array  beamformers.  Furthermore,  the 
proposed  generic  processing  structure  integrates  passive  and  active  modes  of  operation  into 
a  single  signal  processing  scheme. 


3.1  Generic  Multi-Dimensional  Conventional 
Beamforming  Structure 

3.1.1  Line-Array  Conventional  Beamformer 

Consider  an  A-sensor  line  array  receiver  with  uniform  sensor  spacing  8, 
shown  in  Figure  5,  receiving  plane-wave  arrivals  with  direction  of 
propagation  8 .  Then,  as  a  follow  up  of  the  parameter  definition  in  Section  2, 

Tn(0)  =  (n-Y)S  COS0  /  c  ,  (23) 

is  the  time  delay  between  the  1st  and  the  nth  sensor  of  the  line  array  for  an 
incoming  plane  wave  with  direction  9,  as  this  is  illustrated  in  Figures  5  and  6. 

rf„(/„0)  =  ex ,  (24) 
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is  the  n  term  of  the  steering  vector  D  (  f,  0 )  .  Moreover,  because  of 
relations  (16)  and  (17)  the  plane  wave  response  of  the  TV-sensor  line  array 
steered  at  a  direction  0S  can  be  expressed  by, 

B{f,es)  =  D\f,es)x(f).  (25) 

Previous  studies  [1]  have  shown  that  for  a  single  source  this  conventional 
beamformer  without  shading  is  an  optimum  processing  scheme  for  bearing 
estimation.  The  side  lobe  structure  can  be  suppressed  at  the  expense  of  a 
beam  width  increase  by  applying  different  weights  (i.e.,  spatial  shading 
window)  [36].  The  angular  response  of  a  line-array  is  ambiguous  with 
respect  to  the  angle  6s ,  responding  equally  to  targets  at  angle  6 v  and  -6s 
where  6 s  varies  over  [ 0,  ^r] . 


Eq.  (25)  is  basically  a  mathematical  interpretation  of  Figure  6  and  shows  that 
a  line  array  is  basically  a  spatial  filter  because  by  steering  a  beam  in  a 
particular  direction  we  spatially  filter  the  signal  coming  from  that  direction, 
as  this  is  illustrated  in  Figure  5.  On  the  other  hand,  Eq.  (25)  is  fundamentally 
a  discrete  Fourier  transform  relationship  between  the  hydrophone  weightings 
and  the  beam  pattern  of  the  line  array  and  as  such  it  is  computationally  a  very 
efficient  operation.  Flowever,  Eq.  (25)  can  be  generalized  for  non-linear  2- 
dimensional  and  3-dimensional  arrays  and  this  is  discussed  in  the  next 
section. 

As  an  example,  let  us  consider  a  distant  monochromatic  source.  Then  the 
plane  wave  signal  arrival  from  the  direction  6  received  by  a  A -hydrophone 
line  array  is  expressed  by  Eq.  (24).  The  beam  power  pattern  P((6J  is  given 
by  P(f,0s)  -  B{f,6s )  X  B  (/, Q, )  that  takes  the  form 


P(fA)  =  ttxn(f)X:Xf)ex  p 

n= 1  m-\ 


J27tfSnm  COS  ft. 

c 


(26) 


26 


DR  DC  Toronto  TR  2002-058 


where  8}im  is  the  spacing  8(n-m)  between  the  and  m(h  sensors.  As  a 
result  of  Eq.  (26),  the  expression  for  the  power  beam  pattern  P(f  9S),  is 
reduced  to: 


P(JA)  =  \ 


sin 


k8 

N  —  (sin#s  -sin<2) 
X 


sin 


k8 


(sin  9S  -  sin  6) 


(27) 


Let  us  consider  for  simplicity  the  source  bearing  0  to  be  at  array  broadside, 
S=A/2  and  L  =  (N-l)S  is  the  array  size.  Then  Equation  (27)  is  modified  as 
[4,10]: 


P(fA) 


N2  sin2 

kL  sin  9S 

,  A 

'  nL  sin  0S' 

2 

V  A  , 

(28) 


which  is  the  farfield  radiation  or  directivity  pattern  of  the  line  array  as 
opposed  to  near  field  regions.  The  results  in  Equations  (27)  and  (28)  are  for  a 
perfectly  coherent  incident  acoustic  signal  and  an  increase  in  array  size  L 
results  in  additional  power  output  and  a  reduction  in  beamwidth,  which  are 
similar  arguments  with  those  associated  with  the  CRLB  analysis  expressed  by 
Eq.  (21).  The  side-lobe  structure  of  the  directivity  pattern  of  a  line  array, 
which  is  expressed  by  Eq.  (27),  can  be  suppressed  at  the  expense  of  a 
beamwidth  increase  by  applying  different  weights.  The  selection  of  these 
weights  will  act  as  spatial  filter  coefficients  with  optimum  performance 
[5,17,18],  There  are  two  different  approaches  to  select  the  above  weights: 
pattern  optimization  and  gain  optimization.  For  pattern  optimization  the 
desired  array  response  pattern  P(f,  9S)  is  selected  first.  A  desired  pattern  is 
usually  one  with  a  narrow  main  lobe  and  low  sidelobes.  The  weighting  or 
shading  coefficients  in  this  case  are  real  numbers  from  well  known  window 
functions  that  modify  the  array  response  pattern.  Harris'  review  [36]  on  the 
use  of  windows  in  discrete  Fourier  transforms  and  temporal  spectral  analysis 
is  directly  applicable  in  this  case  to  spatial  spectral  analysis  for  towed  line 
array  applications. 

ETsing  the  approximation  sin  9=  9  for  small  9  at  array  broadside,  the  first 
null  in  Eq.  (25)  occurs  at  nLsin9/k=n  or  A9xL/X  =  1.  The  major  conclusion 
drawn  here  for  line  array  applications  is  that  [4,10]: 

AG~X!L  and  AfxT=l  (29) 
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where  T-M s  /Fs  is  the  sensor  time  series  length.  Both  the  above  relations  in 
Eq.  (29)  express  the  well  known  temporal  and  spatial  resolution  limitations  in 
line  array  applications  that  form  the  driving  force  and  motivation  for  adaptive 
and  synthetic  aperture  signal  processing  that  we  will  discuss  later. 

An  additional  constraint  for  sonar  and  ultrasound  applications  requires  that 
the  frequency  resolution  Af  of  the  hydrophone  time  series  for  spatial  spectral 
analysis  that  is  based  on  FFT  beamforming  processing  must  be  such  that 

Af  X  — «1  (30) 

c 

in  order  to  satisfy  frequency  quantization  effects  associated  with  discrete 
frequency  domain  beamforming  following  the  FFT  of  sensor  data  [17,42], 
This  is  because,  in  conventional  beamforming  Finite-duration  Impulse 
Response  (FIR)  filters  are  used  to  provide  realizations  in  designing  digital 
phase  shifters  for  beam  steering.  Since  fast-convolution  signal  processing 
operations  are  part  of  the  processing  flow  of  a  sonar  signal  processor,  the 
effective  beamforming  filter  length  needs  to  be  considered  as  the  overlap  size 
between  successive  snapshots.  In  this  way,  the  overlap  process  will  account 
for  the  wraparound  errors  that  arise  in  the  fast-convolution  processing  [  1 ,40- 
42].  It  has  been  shown  [42]  that  an  approximate  estimate  of  the  effective 
beamforming  filter  length  is  provided  by  Eqs.  (28)  and  (30). 

Because  of  the  linearity  of  the  conventional  beamforming  process,  an  exact 
equivalence  of  the  frequency  domain  narrowband  beamformer  with  that  of 
the  time-domain  beamformer  for  broadband  signals  can  be  derived  [42,43], 
Based  on  the  model  of  Figure  3,  the  time-domain  beamformer  is  simply  a 
time  delaying  [43]  and  summing  process  across  the  hydrophones  of  the  line 
array,  which  is  expressed  by, 

A' 

6(0, .0=  -O  •  (31) 

n  =  1 


Since, 


b(Os,t0  =  l¥FT  {B(fes)}  (32) 

by  using  FFTs  and  fast  convolution  procedures,  continuous  beam-time 
sequences  can  be  obtained  at  the  output  of  the  frequency  domain  beamformer 
[42],  This  is  a  very  useful  operation  when  the  implementation  of 
beamforming  processors  in  sonar  systems  is  considered. 

The  beamforming  operation  in  Eq.  (31)  is  not  restricted  only  for  plane  wave 
signals.  More  specifically,  consider  an  acoustic  source  at  the  near  field  of  a 
line  array  with  rs  the  source  range  and  0  its  bearing.  Then  the  time  delay 
for  steering  at  6  is 
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^  =  (rs  +dL~2  rsdnm  cos  9  )'  2  /  c  • 


(33) 


As  a  result  of  Eq.  (33),  the  steering  vector  dn(f,  9S)  =  exp[j2jtfxs]  will 
include  two  parameters  of  interest,  the  bearing  9  and  range  rs  of  the  source. 
In  this  case  the  beamformer  is  called  focussed  beamformer,  which  is  used 
mainly  in  ultrasound  system  applications.  There  are,  however,  practical 
considerations  restricting  the  application  of  the  focused  beamformer  in 
passive  sonar  line  array  systems  and  these  have  to  do  with  the  fact  that 
effective  range  focussing  by  a  beamformer  requires  extremely  long  arrays. 


Figure  7.  Geometric  configuration  and  coordinate  system  for  a  circular  array  of  sensors. 


3.1.2  Circular  Array  Conventional  Beamformer 

Consider  M- sensors  distributed  uniformly  on  a  ring  of  radius  R  receiving 
plane-wave  arrivals  at  an  azimuth  angle  0  and  an  elevation  angle  (j)  as 
shown  in  Figure  7.  The  plane-wave  response  of  this  circular  array  for  azimuth 
steering  9S  and  an  elevation  steering  (j)s  can  be  written  as  follows: 

=  D\f,GMW{6s)X{f)  ,  (34) 


where  D{f,Os,(j)s )  is  the  steering  vector  with  the  nth  term  being  expressed 

by  dn  (/,  0S  ,0,)  =  expO'2^7?  sin  <j)s  cos (0S  -  Gn)lc )  and 

th 

9m  =  2.mn  I M  is  the  angular  location  of  the  m  sensor  with 
m  =  0,1, ...M  —  1 .  W(9.)  is  a  diagonal  matrix  with  the  off  diagonal  terms 
being  zero  and  the  diagonal  terms  being  the  weights  of  a  spatial  window  to 


DRDC  Toronto  TR  2002-058 


29 


reduce  the  side-lobe  structure  [36],  This  spatial  window,  in  general,  is  not 
uniform  and  depends  on  the  sensor  location  (8m  )  and  the  beam  steering 
direction  (8 v ).  The  beam  power  pattern  P(f  ,6S ,  (j)s )  is  given  by 
P(f,  0S ,  <ps )  =  B(f.  Gs ,  <j>s  )xB  (/,  9S ,  (f)s ) .  The  azimuth  angular  response 
of  the  circular  array  covers  the  range  [0  ,  2ft]  and  therefore  there  is  no 

ambiguity  with  respect  to  the  azimuth  angle  8 . 


3.2  Multidimensional  (3-D)  Array  Conventional 
Beamformer 

Presented  in  this  section  is  a  generic  approach  to  decompose  the  planar,  cylindrical  and 
spherical  array  beamformers  into  coherent  sub-sets  of  line  and/or  circular  array 
beamformers.  In  this  report,  we  will  restrict  the  discussion  to  3D  arrays  with  cylindrical  and 
planar  geometric  configuration.  The  details  of  the  decomposition  process  for  spherical  arrays 
are  similar  and  can  be  found  in  [7,25-28], 

3.2.1  Decomposition  Process  for  2-D  &  3-D  Sensor  Array 
Beamformers 

3.2. 1 . 1  Cylindrical  Array  Beamformer 

Consider  the  cylindrical  array  shown  in  Fig.  8  with  X  sensors  and 
N  =  NM ,  where  N  is  the  number  of  circular  rings  and  M  is  the 
number  of  sensors  on  each  ring.  The  angular  response  of  this 
cylindrical  array  to  a  steered  direction  at  {8S ,  <px )  can  be 
expressed  as 

B(f,  e,  ,*,)  =  ££  Kmx,m  )  o$> 

r=0  ni~  0 

th 

where  wr  m  is  the  (r,m)  term  of  a  3-D  spatial  window,  Xr  m(f) 

’  th 

is  the  (r,m)  term  of  the  matrix  £(  f) ,  ox  X  r  m  (f)  is  the  Fourier 

th  th 

transform  of  the  signal  received  by  the  m  sensor  on  the  r  ring 
and 

dr,m  (/>  >  <Ps )  =  ox]){j2ftf[(r8:  cos  (ps  +  R  sin  <px  cos (8S  -  8m )  /  c)J 

is  the  (r,m)  steering  term  of  D(  J,  0S,  <j>s ).  R  is  the  radius  of 
the  ring,  8.  is  the  distance  between  each  ring  along  z-axis,  r  is 
the  index  for  the  r  ring  and  8  —  2mn  /  M  , 
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m  =  0,1,...,  M  —  1 .  Assuming  wr,m  -  wr  x  wm,  Eq.  (35)  can  be  re¬ 
arranged  as  follows: 


AM 


r= 0 


M- 1 


5X,  t/>X  (/.»,.  A) 


w=0 


(36) 


where  d,  (/,  0,  ,0,)  =  exp[/'2^(r  cos  ^  /  c|J  is  the  r 

term  of  the  steering  vector  for  line-array  beamforming,  wr  is  the 
th 

r  term  of  a  spatial  window  for  line  array  spatial  shading, 

dm (/> 0,A)  =  exp{/2;gf  (/? sin <ps  cos (0,  - 0m ) / c)}  is  the  m* 

term  of  the  steering  vector  for  a  circular  beamformer,  discussed  in 

th 

Section  3.1,  and  wm  is  the  m  term  of  a  spatial  window  for 
circular  array  shading.  Thus,  Eq.  (36)  suggests  the 
decomposition  of  the  cylindrical  array  beamformer  into  two  steps, 
which  is  a  well-known  process  in  array  theory.  The  first  step  is  to 
perform  circular  array  beamforming  for  each  of  the  N  rings  with 
M  sensors  on  each  ring.  The  second  step  is  to  perform  line  array 
beamforming  along  the  z-axis  on  the  A-beam  time  series  outputs 
of  the  first  step.  This  type  of  implementation,  which  is  based  on 
the  decomposition  of  the  cylindrical  beamformer  into  line  and 
circular  array  beamformers  is  shown  in  Figure  8.  The  coordinate 
system  is  identical  to  that  shown  in  Figure  6.  The  decomposition 
process  of  Eq.  (36)  makes  also  the  design  and  incorporation  of  3- 
D  spatial  windows  much  simpler.  Non-uniform  shading  windows 
can  be  applied  to  each  circular  beamformer  to  improve  the  angular 
response  with  respect  to  the  azimuth  angle,  6 .  A  uniform  shading 
window  can  then  be  applied  to  the  line  array  beamformer  to 
improve  the  angular  response  with  respect  to  the  elevation  angle, 

(j) .  Moreover,  the  decomposition  process,  shown  in  Figure  8, 
leads  to  an  efficient  implementation  in  computing  architectures 
based  on  the  following  two  factors: 


•  The  number  of  sensors  for  each  of  these  circular  and  line 
array  beamformers  is  much  less  than  the  total  number  of 
sensors,  N ,  of  the  cylindrical  array.  The  proposed 
decomposition  process  for  the  3-D  beamformer  eliminates 
the  need  for  very  large  memory  and  CPU's  with  very  high 
throughput  requirements  in  one  board  for  real  time 
system  applications. 

•  All  the  circular  and  line  array  beamformers  can  be 
executed  in  parallel,  which  allows  their  implementations 
in  much  simpler  parallel  architectures  with  simpler  CPU's, 
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which  is  a  practical  requirement  for  real  time  system 
applications. 

Thus,  under  the  restriction  w,  m  -  wr  x  w,„  for  3-D  spatial  shading, 
the  decomposition  process  provides  equivalent  beam  time  series 
as  with  those  that  would  have  been  provided  by  a  3-D  cylindrical 
beamformcr,  as  this  is  shown  by  Eqs.  (35)  &  (36). 


3.2. 1.2  Planar  Array  Beamformer 


Consider  the  discrete  planar  array  shown  in  Figure  4  with  N 
sensors  where  N  =  NM  and  M,  N  are  the  number  of  sensors 
along  x-axis  and  y-axis,  respectively.  The  angular  response  of  this 
planar  array  to  a  steered  direction  ( , (f>s )  can  be  expressed  as, 


A'— 1 A/-1 

V  ' 


BtfAM = YYk  JC,  . 


(37) 


;=0  m= 0 


32 


DR  DC  Toronto  TR  2002-058 


elevation  bearing  estimates)  beamformers. 
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where  wr  m  is  the  (r,m)th  term  of  matrix  W(0,(f>)  including  the 
weights  of  a  2-D  spatial  window,  X,  m  (f)  is  the  (r,m) lh 
term  of  the  matrix  X_{  f)  including  the  Fourier  transform  of  the 
received  signal  by  the  (m,r)th  sensor  along  x-axis  and  y-axis, 
respectively.  £)(/,(?  ,0v)  is  the  steering  matrix  having  its  (r,m)th 
term  defined  by 

dr  m  (/,  0s , <ps )  =  exp(  jlnf  ( mSx  sin  0s  +  rSv  cos  0S  cos  (p 


Assuming  that  the  matrix  of  spatial  shading  (weighting)  W(0,</>)  is 
separable  (i.e.,  =  tF,(0)tK,(0)),  Eq.  37  can  be  simplified  as 

follows: 


(38a) 


where,  dr  (/, 0  ,  <ps )  =  exp{j27tfr8y  cos 6S  cos (j)s  I c) ,  is  the  rlh  term 
of  the  steering  vector,  Dy  (/,  6S ,  (j)y )  and 
dm  (/,  Gs ,  <ps )  =  exp(jl7tfmSx  sin  0S  /  c)  is  the  mlh  term  of  the 
steering  vector,  Dx  (/, 6s ,  (f>s ) .  The  summation  term  enclosed  by 

parenthesis  in  Eq.  (38a)  is  equivalent  to  the  response  of  a  line 
array  beamformer  along  x-axis  .  Then  all  the  steered  beams  from 
this  summation  term  form  a  vector  denoted  by  By(f,ds)  .  This 

vector  defines  a  line  array  with  directional  sensors,  which  are  the 
beams  defined  by  the  second  summation  process  of  Eq.  (38a). 
Therefore  Eq.  (38a)  can  be  expressed  as: 

B{.f,ds,(Ps )  =  D'M.es^s  )W^6)Bv(f\9j.  (38b) 

Eq.  (38b)  suggests  that  the  2-D  planar  array  beamformer  can  be 
decomposed  into  two  line  array  beamforming  steps.  The  first  step 
includes  a  line-array  beamforming  along  x-axis  and  will  be 

repeated  N-  times  to  get  the  vector  By(f,Os)  that  includes  the 
beam  times  series  b,(f  9J  ,  where  the  index  r  —  0,1,...,  N  —  1  is 
along  the  axis  y.  The  second  step  includes  line  array  beamforming 
along  y-axis  and  will  be  done  only  once  by  treating  the  vector 
By(f,0s)  as  the  input  signal  for  the  line  array  beamformer  to 

get  the  output  B{f  ,0 ,,<py)  ■  The  separable  spatial  windows  can 

now  be  applied  separately  on  each  line-array  beamformer  to 
suppress  sidelobe  structures.  Figure  4  shows  the  involved  steps  of 
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decomposing  the  2-D  planar  array  beamformer  into  two  steps  of 
line-array  beamformers.  The  line  array  coordinate  system  is 
identical  with  that  shown  in  Figure  6.  The  decomposition  of  the 
planar  array  beamformer  into  these  two  line-array  beamforming 
steps  leads  to  an  efficient  implementation  based  on  the  following 
two  factors.  First,  the  number  of  the  involved  sensors  for  each  of 
these  line  array  beamformers  is  much  less  than  the  total  number  of 
sensors,  X  ,  of  the  planar  array.  The  proposed  decomposition 
process  for  the  2-D  beamformer  eliminates  the  need  for  very  large 
memory  and  CPU's  with  very  high  throughput  requirements  in  one 
board  for  real  time  system  applications.  Secondly,  all  these  line 
array  beamformers  can  be  executed  in  parallel,  which  allows  their 
implementation  in  much  simpler  parallel  architectures  with 
simpler  CPU's,  which  is  a  practical  requirement  for  real  time 
system  application.  Besides  the  advantage  of  the  efficient 
implementation,  the  proposed  decomposition  approach  makes  the 
application  of  the  spatial  window  much  simpler  to  be 
incorporated. 


3.3  Influence  of  the  Medium’s  Propagation 
Characteristics  on  the  Performance  of  a 
Receiving  Array 

In  ocean  acoustics  and  medical  ultrasound  imaging,  the  wave  propagation  problem  is  highly 
complex  due  to  the  spatial  properties  of  the  non-homogeneous  underwater  and  human  body 
mediums.  For  stationary  source  and  receiving  arrays,  the  space  time  properties  of  the 
acoustic  pressure  fields  include  a  limiting  resolution  imposed  by  these  mediums.  This 
limitation  is  due  either  to  the  angular  spread  of  the  incident  energy  about  a  single  arrival  as  a 
result  of  the  scattering  phenomena,  or  to  the  multipaths  and  their  variation  over  the  aperture 
of  the  receiving  array. 

More  specifically,  an  acoustic  signal  that  propagates  through  anisotropic  mediums  will 
interact  with  the  transmitting  medium  microstructure  and  the  rough  boundaries,  resulting  in 
a  net  field  that  is  characterized  by  irregular  spatial  and  temporal  variations.  As  a 
consequence  of  these  interactions,  a  point  source  detected  by  a  high-angular  resolution 
receiver  is  perceived  as  a  source  of  finite  extent.  It  has  been  suggested  [47]  that  due  to  the 
above  spatial  variations  the  sound  field  consists  not  of  parallel,  but  of  superimposed 
wavefronts  of  different  directions  of  propagation.  As  a  result,  coherence  measurements  of 
this  field  by  a  receiving  array  give  an  estimate  for  the  spatial  coherence  function.  In  the 
model  for  the  spatial  uncertainty  of  the  above  study  [47],  the  width  of  the  coherence  function 
is  defined  as  the  coherence  length  of  the  medium  and  its  reciprocal  value  is  a  measure  of  the 
angular  uncertainty  caused  by  the  scattered  field  of  the  underwater  environment. 

By  the  coherence  of  acoustic  signals  in  the  sea  or  in  the  human  body,  we  mean  the  degree  to 
which  the  acoustic  pressures  are  the  same  at  two  points  in  the  medium  of  interest  located  a 
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given  distance  and  direction  apart.  Pressure  sensors  placed  at  these  two  points  will  have 
phase  coherent  outputs  if  the  received  acoustic  signals  are  perfectly  coherent;  if  the  two 
sensor  outputs,  as  a  function  of  space  or  time,  are  totally  dissimilar,  the  signals  are  said  to  be 
incoherent.  Thus,  the  loss  of  spatial  coherence  results  in  an  upper  limit  on  the  useful  aperture 
of  a  receiving  array  of  sensors  [10].  Consequently,  knowledge  of  the  angular  uncertainty  of 
the  signal  caused  by  the  medium  is  considered  essential  in  order  to  determine  quantitatively 
the  influence  of  the  medium  on  the  array  gain,  which  is  also  influenced  significantly  by  a 
partially  directive  anisotropic  noise  background.  Therefore,  for  a  given  non-isotropic 
medium,  it  is  desirable  to  estimate  the  optimum  array  size  and  achievable  array  gain  for 
sonar  and  ultrasound  array  applications. 

For  geometrical  simplicity  and  without  any  loss  of  generality  we  consider  the  case  of  a 
receiving  line  array.  Quantitative  estimates  of  the  spatial  coherence  for  a  receiving  line 
array  are  provided  by  the  cross  spectral  density  matrix  in  frequency  domain  between  any  set 
of  two  sensor  time  series  of  the  line  array.  An  estimate  of  the  cross  spectral  density  matrix 
R(f)  with  its  nmth  term  defined  by 

R,(f,SJ  =  E[xx,r)XHf)]  '  (39) 


The  above  space-frequency  correlation  function  can  be  related  to  the  angular  power 
directivity  pattern  of  the  source,  %(f, 6 ),  via  a  Fourier  transformation  by  using  a 

generalization  of  Bello's  concept  [48]  of  time-frequency  correlation  function  [/  <=>  2 ivf  ] 
into  space  [<5m  <44>  2;z/’sin#/c]  ,  which  gives 


nil 


K,XfAJ=  Jy,(/,0)exp 

-n!  2 


[-j2nfS,jn 
c  J 


dd 


(40) 


or 


s's;  2 


^ ,  (/,»)=  J  K.  (/.  S„„ )  exp02^"e]rffc. ) 


-A’  5/2 


J 


(41) 


The  above  transformation  can  be  converted  into  the  following  summation: 


C/2 

/U/,,O=A0  I  >r,(/„A)exp 


g=-G'2 


r~/2  nfJXm  sin(gA(?)l 

L  c  J 


cos(gA  6) 


(42) 


where  A9  is  the  angle  increment  for  sampling  the  angular  power  directivity  pattern,  9g=gA9, 
g  is  the  index  for  the  samples  and  G  is  the  total  number  of  samples. 

For  line  array  applications,  the  power  directivity  pattern  (calculated  for  a  homogeneous  free 
space)  due  to  a  distant  source,  which  is  treated  as  a  point  source,  should  be  a  delta  function. 
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Estimates,  however,  of  the  source's  directivity  from  a  line  array  operating  in  an  anisotropic 
ocean  are  distorted  by  the  underwater  medium.  In  other  words,  the  directivity  pattern  of  the 
received  signal  is  the  convolution  of  the  original  pattern  and  the  angular  directivity  of  the 
medium  (i.e.,  the  angular  scattering  function  of  the  underwater  environment).  As  a  result  of 
the  above,  the  angular  pattern  of  the  received  signal,  by  a  receiving  line  array  system,  is  the 
scattering  function  of  the  medium. 

In  this  report,  the  concept  of  spatial  coherence  is  used  to  determine  the  statistical  response  of 
a  line  array  to  the  acoustic  field.  This  response  is  the  result  of  the  multipath  and  scattering 
phenomena  discussed  before,  and  there  are  models  [10,47]  to  relate  the  spatial  coherence 
with  the  physical  parameters  of  an  anisotropic  medium  for  measurement  interpretation.  In 
these  models,  the  interaction  of  the  acoustic  signal  with  the  transmitting  medium  is 
considered  to  result  in  superimposed  wavefronts  of  different  directions  of  propagation. 

Then  Eqs.  (24),  (25),  which  define  a  received  sensor  signal  from  a  distant  source,  are 
expressed  by 

x„  (0  =  Tj a,  exp[-/ 2  iff,  (t,  -  -  9, )  j  +  eni  (0,  o e )  (43) 


where  ,  and./  is  the  number  of  superimposed  waves.  As  a  result,  a 

generalized  form  of  the  crosscorrelation  function  between  two  sensors,  which  has  been 
discussed  by  Carey  and  Moseley  [10],  is 


X2(f)exp 


r  r*  vi 


V  LJ 


k-\,  or  1.5  or  2, 


(44) 


where  Lc  is  the  correlation  length  and  X2  (f)  is  the  mean  acoustic  intensity  of  a  received 
sensor  time  sequence  at  the  frequency  bin  f.  A  more  explicit  expression  for  the  Gaussian 
form  of  Eq.  (44)  is  given  in  [47], 


K.UA.)- ?(/)<* 

and  the  crosscorrelation  coefficients  are  given  from 

At  the  distance  Lc  =  c  /  (2  7lf(7h ),  called  "the  coherence  length  ",  the  correlation  function 
in  Eq.  (46)  will  be  0.6.  This  critical  length  is  determined  from  experimental  coherence 
measurements  plotted  as  a  function  of  8nm.  Then  a  connection  between  the  medium's  angular 
uncertainty  and  the  measured  coherence  length  is  derived  as 


/  2 


(45) 


(46) 
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=  i  /  4. 


and 


L  =2nd]nuf  /  c 


(47) 


here  8\m  is  the  critical  distance  between  the  first  and  the  m  th  sensors  at  which  the 
coherence  measurements  get  smaller  than  0.6.  Using  the  above  parameter  definition,  the 
effective  aperture  size  and  array  gain  of  a  deployed  towed  line  array  can  be  determined 
[10,47]  for  a  specific  underwater  ocean  environment. 

Since  the  correlation  function  for  a  Gaussian  acoustic  field  is  given  by  Eq.  (45),  the  angular 
scattering  function  &(f,9)  of  the  medium  can  be  derived.  Using  Eq.  (41)  and  following  a 
rather  simple  analytical  integral  evaluation,  we  have 


®(f..0)  = 


(48) 


where  oq~  c/(27ifSlini).  This  is  an  expression  for  the  angular  scattering  function  of  a 
Gaussian  underwater  ocean  acoustic  field  [10,47], 

It  is  apparent  from  the  above  discussion  that  the  estimates  of  the  cross-correlation 
coefficients  pnm  (/] ,  8nm )  are  necessary  in  order  to  define  experimentally  the  coherence 

length  of  an  undewater  or  human  body  medium.  For  details  on  experimental  studies  on 
coherence  estimation  for  undewater  sonar  applications  the  reader  may  review  the  references 
[10,30], 


3.4  Array  Gain 

The  performance  of  an  array  of  sensors  to  an  acoustic  signal  embodied  in  a  noise  field  is 
characterized  by  the  "array  gain  "  parameter,  AG.  The  mathematical  relation  of  this 
parameter  is  defined  by 


XZa„„(.a<d 

AG  -  IOlog-r—f- - 


(49) 


where  pnm  (/] , 8nm )  and  pE>nm(f,8nnJ  denote  the  normalized  crosscorrelation  coefficients  of 

the  signal  and  noise  field,  respectively.  Estimates  of  the  correlation  coefficients  are  given 
from  Eq.  (46).  If  the  noise  field  is  isotropic  that  it  is  not  partially  directive,  then  the 

N  N 

denominator  in  Eq.  (49)  is  equal  to  N,  (i.e.  12 A-JfAJ  =  N  ), 

n=  I  m= 1 

because  the  non  diagonal  terms  of  the  cross-correlation  matrix  for  the  noise  field  are 
negligible.  Then  Equation  (49)  simplifies  to, 
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(50) 


XZ/u/.o 

AG  =  lOlog-2^1-2^ - 

N 

For  perfect  spatial  coherence  across  the  line  array  the  normalized  crosscorrelation 
coefficients  are  pnm(f,  Snm)  =  1  and  the  expected  values  of  the  array  gain  estimates  are, 

AG  =  10  X  log  N .  For  the  general  case  of  isotropic  noise  and  for  frequencies  smaller  than 
the  towed  array's  design  frequency  the  array  gain  term  AG  is  reduced  to  the  quantity  called 
Directivity  Index, 


D/  =  10xlog[(Af-l>y/(A/2)].  (51) 

When  S«  X  and  the  conventional  beamforming  processing  is  employed,  Eq.  (29) 
indicates  that  the  deployment  of  very  large  aperture  arrays  is  required  in  order  to  achieve 
sufficient  array  gain  and  angular  resolution  for  precise  bearing  estimates.  Practical 
deployment  considerations,  however,  usually  limit  the  overall  dimensions  of  an  array  of 
sensors.  In  addition,  the  medium's  spatial  coherence  [10,  30]  sets  an  upper  limit  on  the 
effective  array  length.  In  general,  the  medium's  spatial  coherence  length  is  of  the  order  of 
O(102)?i  [10,30], 

Alternatives  to  large  aperture  arrays  are  signal  processing  schemes  discussed  in  [1]. 
Theoretical  and  experimental  investigations  have  shown  that  bearing  resolution  and 
detectability  of  weak  signals  in  the  presence  of  strong  interferences  can  be  improved  by 
applying  non-conventional  beamformers  such  as  adaptive  beamforming  [1-5,17-24],  or 
acoustic  synthetic  aperture  processing  [1,1 1-16]  to  the  sensor  time  series  of  deployed  short 
ultrasound  arrays,  which  are  discussed  in  the  next  section. 
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4.  Advanced  Beamformers 


4.1  Synthetic  Aperture  Processing 

Various  synthetic  aperture  techniques  have  been  investigated  to  increase  signal  gain  and 
improve  angular  resolution  for  line  array  systems.  There  is  a  difference,  however,  between 
the  synthetic  aperture  concept  of  ultrasound  systems  with  that  for  sonar  and  radar 
applications.  In  particular,  the  concept  of  synthetic  aperture  processing  in  sonar  and  radar 
system  applications  exploits  the  movement  of  an  array  of  sensors  to  synthesize  a  longer 
aperture.  In  ultrasound  system  applications,  the  synthetic  aperture  processing  concept 
integrates  the  acquired  data  from  sub-apertures  of  a  larger  physical  array.  Thus,  for 
ultrasound  systems,  the  synthetic  aperture  concept  is  simpler  since  there  is  no  movement  of 
the  array,  while  for  radar  and  sonar  systems  the  movement  of  the  array  introduces  numerous 
errors  and  artifacts.  In  what  follows,  we  will  review  synthetic  aperture  techniques  for  sonar 
applications  that  have  been  tested  successfully  with  real  data  [1 1-16].  The  same  techniques 
can  be  implemented  successfully  in  ultrasound  systems  and  this  will  be  demonstrated  later  in 
the  section  of  this  report  that  includes  results  from  numerical  simulations  and  experimental 
data. 

Let  us  start  with  the  a  few  theoretical  remarks.  The  plane  wave  response  of  a  line  array  to  a 
distant  monochromatic  signal,  received  by  the  nth  element  of  the  array,  is  expressed  by  Eqs. 
(23),  (24)  and  (25).  In  the  above  expressions,  the  frequency  /  includes  the  Doppler  shift 
due  to  a  combined  movement  of  the  receiving  array  and  the  source  (or  object  reflecting  the 
incoming  acoustic  wavefront)  radiating  signal.  Let  v,  denote  the  relative  speed;  it  is 
assumed  here  that  the  component  of  the  source’s  velocity  along  its  bearing  is  negligible.  If 
f  ,  is  the  frequency  of  the  stationary  field,  then  the  frequency  of  the  received  signal  is 
expressed  by 


/  =  /0(1±ysin#/c) 


(52) 


and  an  approximate  expression  for  the  received  sensor  time  series  (18)  and  (43)  is  given  by 


x„  (/,  )=  A  exp 


j  2tf<, 


.  Vti+(n- 1)£  . 

t - - - sin  6 

c  ) 


+  £., 


(53) 


r  seconds  later,  the  relative  movement  between  the  receiving  array  and  the  radiated  source  is 
VT .  By  proper  choice  of  the  parameters  V  and  r,  we  have  vv=qS,  where  q  represents  the 
number  of  sensor  positions  that  the  array  has  moved,  and  the  received  signal,  x„(t,+z)  is 
expressed  by, 


c„  (/,.  +  t)=  exp(j2jtf„T)A  exp 


JW„ 


Vt:  +  (q  +  n  — 1)<5  . 


sin  6 


+  £, 


(54) 
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As  a  result,  we  have  the  Fourier  transform  of  x„(ii+r),  as 


X.(f),=e*p(jWs)X.(f) 


(55) 


where,  X n  (/ )  and  X n  (/ )  are  the  DFTs  of  x„(tj+T),  and  xn(t),  respectively.  If  the  phase 

term  exp  (—  j2ltf0'l)  is  used  to  correct  the  line  array  measurements  shown  in  (55),  then 

the  spatial  information  included  in  the  successive  measurements  at  t=t ,  and  /=/,+ris 
equivalent  to  that  derived  from  a  line  array  of  (q+N)  sensors.  When  idealized  conditions  are 
assumed,  the  phase  correction  factor  for  (52)  in  order  to  form  a  synthetic  aperture,  is 
exp(-  j2nfot)  ■  However,  this  phase  correction  estimate  requires  a  priori  knowledge  of 

the  source  receiver  relative  speed,  v  and  accurate  estimates  for  the  frequency  / of  the 
received  signal.  An  additional  restriction  is  that  the  synthetic  aperture  processing 
techniques  have  to  compensate  for  the  disturbed  paths  of  the  receiving  array  during  the 
integration  period  that  the  synthetic  aperture  is  formed.  Moreover,  the  temporal  coherence 
of  the  source  signal  should  be  greater  or  at  least  equal  to  the  integration  time  of  the  synthetic 
aperture. 

1 .  At  this  point  it  is  important  to  review  a  few  fundamental  physical  arguments 
associated  with  passive  synthetic  aperture  processing.  In  the  past  [13]  there  was  a 
conventional  wisdom  regarding  synthetic  aperture  techniques,  which  held  that 
practical  limitations  prevent  them  from  being  applicable  to  real-world  systems.  The 
issues  were  threshold. 

2.  Since  passive  synthetic  aperture  can  be  viewed  as  a  scheme  that  converts  temporal 
gain  to  spatial  gain,  most  signals  of  interest  do  not  have  sufficient  temporal 
coherence  to  allow  a  long  spatially  coherent  aperture  to  be  synthesized. 

3.  Since  past  algorithms  required  a  priori  knowledge  of  the  source  frequency  in  order 
to  compute  the  phase  correction  factor,  as  shown  by  (52)-(55),  the  method  was 
essentially  useless  in  any  bearing  estimation  problem  since  Doppler  would  introduce 
an  unknown  bias  on  the  frequency  observed  at  the  receiver. 

Since  synthetic  aperture  processing  essentially  converts  temporal  gain  to  spatial  gain,,  there 
was  no  “new”  gain  to  be  achieved,  and  therefore,  no  point  to  the  method. 

Recent  work  [12-16]  has  shown  that  there  can  be  realistic  conditions  under  which  all  of 
these  objections  are  either  not  relevant  or  do  not  constitute  serious  impediments  to  practical 
applications  of  synthetic  aperture  processing  in  operational  systems  [1],  Theoretical 
discussions  have  shown  [13]  that  the  above  three  arguments  are  valid  for  cases  that  include 
the  formation  of  synthetic  aperture  in  mediums  with  isotropic  noise  characteristic.  However, 
when  the  noise  characteristics  of  the  received  signal  are  non-isotropic  and  the  receiving 
array  includes  more  than  one  sensor,  then  there  is  spatial  gain  available  from  passive 
synthetic  aperture  processing  and  this  has  been  discussed  analytically  in  [13].  Recently, 
there  have  been  only  two  passive  synthetic  aperture  techniques  [1 1-16,  54]  and  an  MLE 
estimator  [12]  published  in  the  open  literature  that  deal  successfully  with  the  above 
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restrictions.  For  more  details  about  these  techniques  the  reader  may  review  the  references 
[11-16,54,55], 

For  ultrasound  applications,  the  method  that  does  not  require  adaptations  to  address  the 
depth  focusing  requirements  of  ultrasound  systems  is  the  spatial  overlap  correlator  method 
(ETAM  Algorithm),  which  is  discussed  in  detail  in  [1,  11-16]. 

4.2  Adaptive  Beamformers 

Despite  the  geometric  differences  between  the  line  and  circular  arrays,  the  underlying 
beamforming  processes  for  these  arrays,  as  expressed  by  Eqs.  (1 1)  &  (12)  respectively,  are 
time  delay  beamforming  estimators,  which  are  basically  spatial  filters.  However,  optimum 
beamforming  requires  the  beamforming  filter  coefficients  to  be  chosen  based  on  the 
covariance  matrix  of  the  received  data  by  the  N  -sensor  array  in  order  to  optimize  the  array 
response  [15,16],  as  discussed  in  Section  2.  The  family  of  algorithms  for  optimum 
beamforming  that  use  the  characteristics  of  the  noise,  are  called  Adaptive  Beamformers 
[3,17,18,19,20,22,23],  In  this  section  we  will  address  implementation  issues  for  various 
partially  adaptive  variants  of  the  MVDR  method  and  a  GSC  adaptive  beamformer  [1,37], 

Furthermore,  the  implementation  of  adaptive  schemes  in  real  time  systems  is  not  restricted 
into  one  method,  such  as  the  MVDR  technique  that  is  discussed  next.  In  fact,  the  generic 
concept  of  the  sup-aperture  multi-dimensional  array  introduced  in  the  paper  allows  for  the 
implementation  of  a  wide  variety  of  adaptive  schemes  in  operational  systems  [7,25-28].  As 
for  the  implementation  of  adaptive  processing  schemes  in  active  systems,  the  following 
issues  need  to  be  addressed. 

For  applications  that  require  Doppler  processing,  the  outputs  of  the  adaptive  algorithms  are 
required  to  provide  coherent  beam  time  series  to  facilitate  the  post-processing.  This  means 
that  these  algorithms  should  exhibit  near-instantaneous  convergence  and  provide  continuous 
beam  time  series  that  have  sufficient  temporal  coherence  to  correlate  with  the  reference 
signal  in  matched  filter  processing  [1], 

In  a  previous  study  [1],  possible  improvement  in  convergence  periods  of  two  algorithms  in 
the  sub-aperture  configuration  was  investigated.  The  Griffiths-Jim  Generalized  Sidelobe 
Canceller  (GSC)  [18,44]  coupled  with  the  Normalized  Least  Mean  Square  (NLMS)  adaptive 
filter  [45]  has  been  shown  to  provide  near-instantaneous  convergence  under  certain 
conditions  [1,37].  The  GSC/NLMS  in  the  sub-aperture  configuration  w'as  tested  under  a 
variety  of  conditions  to  determine  if  it  could  yield  performance  advantages,  and  if  its 
convergence  properties  could  be  exploited  over  a  wider  range  of  conditions  [1,37].  The 
Steered  Minimum  Variance  Beamformer  (STMV)  is  a  variant  of  the  Minimum  Variance 
Distortionless  Response  (MVDR)  beamformer  [38],  By  applying  namwband  adaptive 
processing  on  bands  of  frequencies,  extra  degrees  of  freedom  are  introduced.  The  number  of 
degrees  of  freedom  is  equal  to  the  number  of  frequency  bins  in  the  processed  band.  In  other 
words,  increasing  the  number  of  frequency  bins  processed  decreases  the  convergence  time 
by  a  corresponding  factor.  This  is  due  to  the  fact  that  convergence  now  depends  on  the 
observation  time  bandwddth  product,  as  opposed  to  observation  time  in  the  MVDR 
algorithm  [38,39], 
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The  STMV  beamformer  in  its  original  form  was  a  broadband  processor.  In  order  to  satisfy 
the  requirements  for  matched  filter  processing,  it  was  modified  to  produce  coherent  beam 
time  series  [1],  The  ability  of  the  STMV  narrowband  beamformer  to  produce  coherent  beam 
time  series  has  been  investigated  in  another  study  [37].  Also,  the  STMV  narrowband 
processor  was  implemented  in  the  sub-aperture  configuration  to  produce  near-instantaneous 
convergence  and  to  reduce  the  computational  complexity  required.  The  convergence 
properties  of  both  the  full  aperture  and  sub-aperture  implementations  have  been  investigated 
for  line  arrays  of  sensors  [1,37]. 


4.2.1  Minimum  Variance  Distortionless  Response  (MVDR) 

The  goal  is  to  optimize  the  beamformer  response  so  that  the  output  contains 
minimal  contributions  due  to  noise  and  signals  arriving  from  directions  other 
than  the  desired  signal  direction.  For  this  optimization  procedure  it  is  desired 
to  find  a  linear  filter  vector  W(fi,&)  which  is  a  solution  to  the  constrained 
minimization  problem  that  allows  signals  from  the  look  direction  to  pass  with 
a  specified  gain  [17,18], 

Minimize:  &MV  =  W*  @)  ,  subject  to 

w*(fi,ejD(fi,0)= i  (56) 

where  JXf^O)  is  the  conventional  steering  vector  based  on  Eq.  (24).  The 
solution  is  given  by. 


(57) 


The  above  solution  provides  the  adaptive  steering  vectors  for  beamforming 
the  received  signals  by  the  N  -hydrophone  line  array.  Then  in  frequency 
domain,  an  adaptive  beam  at  a  steering  9S  is  defined  by 

(58) 

and  the  corresponding  conventional  beams  are  provided  by  Equation  (25). 


4.2.2  Generalized  Sidelobe  Canceller  (GSC) 

The  Generalized  Sidelobe  Canceller  (GSC)  [44]  is  an  alternative  approach  to 
the  MVDR  method.  It  reduces  the  adaptive  problem  to  an  unconstrained 
minimization  process.  The  GSC  formulation  produces  a  much  less 
computationally  intensive  implementation.  In  general  GSC  implementations 
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have  complexity  0(N2),  as  compared  to  0(N1)  for  MVDR  implementations, 
where  N  is  the  number  of  sensors  used  in  the  processing.  The  basis  of  the 
reformulation  of  the  problem  is  the  decomposition  of  the  adaptive  filter 
vector  W(fnO)  into  two  orthogonal  components,  w  and  -v,  where  w  and  v 
lie  in  the  range  and  the  null  space  of  the  constraint  of  Eq.  (56),  such  that 

A  matrix  C  which  is  called  signal  blocking 

matrix,  may  be  computed  from  C  I  =  0  where  /  is  a  vector  of  ones.  This 
matrix  C  whose  columns  form  a  basis  for  the  null  space  of  the  constraint  of 
Eq.  (56)  will  satisfy  y  =  Qj}  ,  where  U  is  defined  below  by  Eq. 
(60).  The  adaptive  filter  vector  may  now  be  defined  as  W  r=w—Gu  and 
yields  the  realization  shown  in  Figure  9.  Then  the  problem  is  reduced  to: 

Minimize:  a]  =  {[w  -  Cu]' R[W  -  Cu])  (59) 

which  is  satisfied  by: 

",  ‘{C'Rcfc’Rw 


Figure  9.  Basic  processing  structure  for  the  memoryless  GSC. 


The  Griffiths-Jim  Generalized  Sidelobe  Canceller  (GSC)  in  combination  with 
the  Normalized  Least  Mean  Square  (NLMS)  adaptive  algorithm  has  been 
shown  to  yield  near  instantaneous  convergence  [44,45].  Figure  9  shows  the 
basic  structure  of  the  so  called  Memoryless  Generalized  Sidelobe  Canceller. 
The  time  delayed  by  Tr:  (  0S )  sensor  time  series  defined  by  Equations  (7), 
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(23)  form  the  pre-steered  sensor  time  series,  which  are  denoted  by 
x n  (ti ,  T„  (0 s )).  In  frequency  domain  these  pre-steered  sensor  data  are 

denoted  by  Xn  (ft ,  0S  )  and  form  the  input  data  vector  for  the  adaptive 
scheme  in  Figure  8.  On  the  left  hand  side  branch  of  this  figure  the 
intermediate  vector  Z  ,  6s  )  is  the  result  of  the  signal  blocking  matrix  C 

being  applied  to  the  input  X  (/] ,  6S ) . 

Next,  the  vector  Z  {ji ,  0s )  is  an  input  to  the  NLMS  adaptive  filter.  The 

output  of  the  right  hand  branch  is  simply  the  shaded  conventional  output. 
Then,  the  output  of  this  processing  scheme  is  the  difference  between  the 
adaptive  filter  output,  and  the  “conventional”  output: 

e(f,  es )  =  b(jt  ,*,)-«*  (/, ,  es  jz<jt ,  es )  (61) 


The  adaptive  filter,  at  convergence,  reflects  the  sidelobe  structure  of  any 
interferers  present,  and  it  is  removed  from  the  conventional  beamformer 
output.  In  the  case  of  the  Normalized  LMS  (NLMS)  this  adaptation  process 
can  be  represented  by: 


Vfc+l 


(fiA)  =  uk(fiA)+—% 


jiXeJfA) 


a+r  UAWA) 


z(fM 


(62) 


where,  k  is  the  iteration  number,  a  is  a  small  positive  number  designed  to 
maintain  stability.  The  parameter  jd  is  the  convergence  controlling  parameter 
or  “step  size”  for  the  NLMS  algorithm. 


4.2.3  Steered  Minimum  Variance  Broadband  Adaptive  (STMV) 

Krolik  and  Swingler  [38]  have  shown  that  the  convergence  time  for  broad¬ 
band  source  location  can  be  reduced  by  using  the  space-time  statistic  called 
the  steered  covariance  matrix  (STCM).  This  method  achieves  significantly 
shorter  convergence  times  than  adaptive  algorithms  that  are  based  on  the 
narrowband  cross  spectral  density  matrix  (CSDM)  [  1 7, 1 8]  without  sacrificing 
spatial  resolution.  In  fact,  the  number  of  statistical  degrees  of  freedom 
available  to  estimate  the  STCM  is  approximately  the  time-bandwidth  product 
( TxBW)  as  opposed  to  the  observation  time,  ( T-MJFS ,  Fs  being  the 
sampling  frequency)  in  CSDM  methods.  This  provides  an  improvement  of 
approximately  B  W,  the  size  of  the  broad-band  source  bandwidth,  in 
convergence  time.  The  conventional  beamformer's  output  in  frequency 
domain  is  shown  by  Eq.  (25).  The  corresponding  time  domain  conventional 
beamformer  output  b(th  6S)  is  the  weighted  sum  of  the  steered  sensor  outputs, 
as  expressed  by  Eq.  (32).  Then,  the  expected  broadband  beam  power  ,B(6)  is 
given  by: 
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(63) 


s(e,  )=E\b(e„  (,)|}=  *•£{?(/,., r,(9))?(/„i-.(e))> 

where  the  vector  includes  the  weights  for  spatial  shading  [36]. 

The  term  6s  )  =  £{*(/, ,  Tn  {0  s  ))x  *  (/, ,  Tm  (0  s  ))]  (64) 


is  defined  as  the  steered  covariance  matrix  (STCM)  in  time  domain  and  is 
assumed  to  be  independent  of  /,  in  stationary  conditions.  The  name  STCM  is 
derived  from  the  fact  that  the  matrix  is  computed  by  taking  the  covariance  of 
the  presteered  time  domain  sensor  outputs.  Suppose X„(fi)  is  the  Fourier 
transform  of  the  sensor  outputs  x„(t)  and  assuming  that  the  sensor  outputs  are 
approximately  band  limited.  Under  these  conditions  the  vector  of  steered  (or 
time  delayed)  sensor  outputs  x„(t„  t„(9J)  can  be  expressed  by 

/+// 

4', .  ))  =  I  rt  </, ,  e,  )x(fk  )exp(/2#p, ) 

*=/  (65) 


where  T(fk,6)  is  the  diagonal  steering  matrix  in  Eq.  (66)  below  with  elements 
identical  to  the  elements  of  the  conventional  steering  vector,  £)(  fj,0) 


1  .  .  .  0 

0  </,(/*  ,0) 


T(.fk,6)  = 


(66) 


0  ...  dK  (/*  ,0) 


Then  it  follows  directly  from  the  above  equations  that 

/+// 

i 'WAhTj(fl.e,)R(f„r(o,) 

*=/  (67) 

where  the  index  k  -1,1+1  ,.,.,1+H  refers  to  the  frequency  bins  in  a  band  of 
interest  4/i  and  R.{fk )  is  the  Cross  Spectral  Density  Matrix  (CSDM)  for  the 

frequency  bin^.  This  suggests  that  <P(Af,9s)  in  Eq.  (64)  can  be  estimated 
from  the  CSDM,  R(f/J  and  T(fk.6)  expressed  by  Eq.  (66).  In  the  steered 
minimum  variance  method  (STMV),  the  broadband  spatial  power  spectral 
estimate  5(7?,)  is  given  by  [38]: 

5(e,)=[r<j>(A/,s,),/['  r68) 
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The  Steered  Minimum  Variance  Algorithm  differs  from  the  basic  MVDR 
algorithm  in  that  the  STMV  algorithm  yields  a  STCM  that  is  composed  from 
a  band  of  frequencies  and  the  MVDR  algorithm  uses  a  CSDM  that  is  derived 
from  a  single  frequency  bin.  Thus,  the  additional  degrees  of  freedom  of 
STMV  compared  to  those  of  CSDM  provide  a  more  robust  adaptive  process. 


However,  estimates  of  B(6)  according  to  Eq.  (68)  do  not  provide  coherent 
beam  time  series,  since  they  represent  the  broadband  beam  power  output  of 
an  adaptive  process.  In  this  investigation  [1]  we  have  modified  the  estimation 
process  of  the  STMV  matrix  in  order  to  get  the  complex  coefficients  of 
0(Af  9S)  for  all  the  frequency  bins  in  the  band  of  interest. 


The  STMV  algorithm  may  be  used  in  its  original  form  to  generate  an  estimate 
of  0(Af,  0)  for  all  the  frequency  bands  Af,  across  the  band  of  the  received 
signal.  Assuming  stationarity  across  the  frequency  bins  of  a  band  Af,  then  the 
estimate  of  the  STMV  may  be  consider  to  be  approximately  the  same  with  the 
narrowband  estimate  <P(f„,  0)  for  the  center  frequency  f,  of  the  band  Af.  In 
this  case,  the  narrowband  adaptive  coefficients  may  be  derived  from 


“'(/».  »)= 


Djr.f) 

D '  (/..«)*(/..  A/,  oy'D(f„  e)' 


(69) 


The  phase  variations  of  w  (/0 , 9)  across  the  frequency  bins  i-l,,l+l,...,l+H 
(where  H  is  the  number  of  bins  in  the  band  Af  ),  are  modeled  by, 

w„(f ,  6)=exp[27fl  W(Af  6)],  i=l,l+l,...,l+H  (70) 

where,  *F„  (Afd)  is  a  time  delay  term  derived  from, 

'Rn(Af,0)  =  F[wn(Af,0),27tfo].  (71) 


Then  by  using  the  adaptive  steering  weights  wn  (Af  6),  that  are  provided  by 
Eq.  (70),  the  adaptive  beams  are  formed  as  shown  by  Eq.  (58).  Figure  10 
shows  the  realization  of  the  STMV  beamformer  and  provides  a  schematic 
representation  of  the  basic  processing  steps  that  include: 

1 .  time  series  segmentation,  overlap  and  FFT.  shown  by  the  group  of 
blocks  at  the  top-left  part  of  the  schematic  diagram. 

2.  Formation  of  steered  covariance  matrix.  [Eqs.  (64),  (67)]  shown  by 
the  two  blocks  at  the  bottom  left  hand  side  of  Figure  10. 
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Figure  10.  Realization  of  the  steered  co variance  adaptive  beamformer. 


3.  Inversion  of  covariance  matrix  using  Cholesky  factorization, 
estimation  of  adaptive  steering  vectors  and  formation  of  adaptive  beams  in 
frequency  domain,  presented  by  the  middle  and  bottom  blocks  at  the  right 
hand  side  of  Figure  10;  and  finally 

4.  formation  of  adaptive  beams  in  time  domain  through  IFFT, 
discarding  of  overlap  and  concatenation  of  segments  to  form  continuous 
beam  time  series,  which  is  shown  by  the  top  right  hand  side  block. 

The  various  indexes  in  Figure  10  provide  details  for  the  implementation  of 
the  STMV  processing  flow  in  a  generic  computing  architecture.  The  same 
figure  indicates  that  estimates  of  the  Steered  Covariance  Matrix  (STCM)  is 
based  on  an  exponentially  weighted  time  average  of  the  current  and  previous 
STCM,  which  is  discussed  in  the  next  section. 
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5.  Implementation  Considerations 

The  conventional  and  adaptive  steering  vectors  for  steering  angles  9 s ,  <f)s  discussed  in 

Section  3  are  integrated  in  a  frequency  domain  beamforming  scheme,  which  is  expressed  by 
Eqs.  (25), (28),  (34)  and  (58).  The  beam  time  series  are  formed  by  Eq.  (32).  Thus,  the 
frequency  domain  adaptive  and  conventional  outputs  are  made  equivalent  to  the  fast  Fourier 
transform  (FFT)  of  the  time  domain  beamforming  outputs  with  proper  selection  of 
beamforming  weights  and  careful  data  partitioning.  This  equivalence  corresponds  to 
implementing  FIR  filters  via  circular  convolution  [40-42]. 

Matrix  inversion  is  another  major  implementation  issue  for  the  adaptive  schemes  discussed 
in  this  report.  Standard  numerical  methods  for  solving  systems  of  linear  equations  can  be 
applied  to  solve  for  the  adaptive  weights.  The  range  of  possible  algorithms  includes: 

Cholesky  factorization  of  the  covariance  matrix  R(fj),  [17,29].  This  allows  the 
linear  system  to  be  solved  by  backsubstitution  in  terms  of  the  received  data 
vector.  Note  that  there  is  no  requirement  to  estimate  the  sample  covariance  matrix 
and  that  there  is  a  continuous  updating  of  an  existing  Cholesky  factorization. 

QR  decomposition  [19,21,24]  of  the  received  vector  X(ft)  ,  that  includes  the 
conversion  of  a  matrix  to  upper  triangular  form  via  rotations.  The  QR 
decomposition  method  has  better  stability  than  the  Cholesky  factorization 
algorithm,  but  it  requires  twice  as  much  computational  efforts  than  the  Cholesky 
approach. 

SVD  (Singular  Value  Decomposition)  method  [  1 9,2 1 ,24],.  This  is  the  most  stable 
factorization  technique.  It  requires,  however,  three  times  more  computational 
requirements  than  the  QR  decomposition  method. 

In  this  implementation  study  we  have  applied  the  Cholesky  factorization  and  the  QR 
decomposition  techniques  in  order  to  get  solutions  for  the  adaptive  weights.  Our  experience 
suggests  that  there  are  no  noticable  differences  in  performance  between  the  above  two 
methods  [1], 

The  main  consideration,  however,  for  implementing  adaptive  schemes  in  real  time  systems 
are  associated  with  the  requirements  derived  from  Eqs.  (57),  (68),  which  require  knowledge 
of  second  order  statistics  for  the  noise  field.  Although  these  statistics  are  usually  not  known, 
they  can  be  estimated  from  the  received  data  [17,18,23]  by  averaging  a  large  number  of 
independent  samples  of  the  covariance  matrixes  R  (f\)  or  by  allowing  the  iteration  process 
of  the  adaptive  GSC  schemes  to  converge  [1,37].  Thus,  if  K  is  the  effective  number  of 
statistically  independent  samples  of  R  (f ),  then  the  variance  on  the  adaptive  beam  output 
power  estimator  detection  statistic  is  inversely  proportional  to  ( K-N+l ),  [17,18,22],  where 
N  is  the  number  of  sensors.  Theoretical  suggestions  [23]  and  our  empirical  observations 
suggest  that  K  needs  to  be  three  to  four  times  the  size  of  N  in  order  to  get  coherent  beam 
time  series  at  the  output  of  the  above  adaptive  schemes.  In  other  words,  for  arrays  with  a 
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large  number  of  sensors,  the  implementation  of  adaptive  schemes  as  statistically  optimum 
beamformers  would  require  the  averaging  of  a  very  large  number  of  independent  samples  of 
R(fi)  in  order  to  derive  an  unbiased  estimate  of  the  adaptive  weights  [23],  In  practice  this 
is  the  most  serious  problem  associated  with  the  implementation  of  adaptive  beamformers  in 
real  time  systems. 

Owsley  [17,29]  has  addressed  this  problem  with  two  important  contributions.  His  first 
contribution  is  associated  with  the  estimation  procedure  of  R  (fi).  His  argument  is  that  in 
practice,  the  covariance  matrix  cannot  be  estimated  exactly  by  time  averaging  because  the 
received  signal  vector  X (/)  is  never  truly  stationary  and/or  ergodic.  As  a  result,  the 
available  averaging  time  is  limited.  Accordingly,  one  approach  to  the  time-varying  adaptive 
estimation  of  R(fi)  at  time  4  is  to  compute  the  exponentially  time  averaged  estimator 
(geometric  forgetting  algorithm)  at  time  4 : 

(/,)+( 1  -Mmf.X  a,)  (72) 

where  /r  is  a  smoothing  factor  ( 0</i<1  )  that  implements  the  exponentially  weighted  time 
averaging  operation.  The  same  principle  has  also  been  applied  in  the  GSC  scheme  [1,37], 
Use  of  this  type  of  exponential  window  to  update  the  covariance  matrix  is  a  very  important 
factor  in  the  implementation  of  adaptive  algorithms  in  real  time  systems. 

Owsley's  [29]  second  contribution  deals  with  the  dynamics  of  the  data  statistics  during  the 
convergence  period  of  the  adaptation  process.  As  mentioned  above,  the  implementation  of 
an  adaptive  beamformer  with  a  large  number  of  adaptive  weights  in  a  large  array  sonar 
system,  requires  very  long  convergence  periods  that  will  eliminate  the  dynamical 
characteristics  of  the  adaptive  beamformer  to  detect  the  time  varying  characteristics  of  a 
received  signal  of  interest.  A  natural  way  to  avoid  temporal  stationarity  limitations  is  to 
reduce  the  number  of  adaptive  weights  requirements.  Owsley's  [29]  sub-aperture 
configuration  for  line  array  adaptive  beamforming  reduces  significantly  the  number  of 
degrees  of  freedom  of  an  adaptation  process.  His  concept  has  been  applied  to  line  arrays,  as 
discussed  in  References  [1,37].  However,  extension  of  the  sub-aperture  line  array  concept 
for  multi-dimensional  arrays  is  not  a  trivial  task.  In  the  following  sections,  the  sub-aperture 
concept  is  generalized  for  circular,  cylindrical,  planar  and  spherical  arrays. 


5.1  Signal  Cancellation  Effects  of  the  Adaptive 
Algorithms 

Testing  of  the  adaptive  algorithms  of  this  study  for  signal  cancellation  effects  was  carried 
out  with  simulations  that  included  two  narrowband  signals  arriving  from  64  degrees  and  69 
degrees  [1,37],  The  simulations  included  full  aperture  and  sub-aperture  (SA) 
implementation,  which  is  discussed  in  the  next  section.  All  of  the  parameters  of  the  signals 
were  set  to  the  same  values  for  all  the  beamformers,  conventional,  GSC/NLMS,  GSC- 
SA/NLMS,  STMV  and  STMV-SA.  Details  about  the  above  simulated  signal  cancellation 
effects  can  be  found  in  Reference  [1,37],  In  the  narrowband  outputs  of  the  conventional 
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beamformer,  the  signals  appear  at  the  frequency  and  beam  at  which  they  were  expected.  As 
anticipated,  however,  the  sidelobes  are  visible  in  a  number  of  other  beams.  The  gram  outputs 
[1,37]  of  the  GSC/STMV  algorithm  indicated  that  there  is  signal  cancellation. 

In  each  case,  the  algorithm  failed  to  detect  either  of  the  two  narrowband  signals.  This 
suggests  that  there  is  a  shortcoming  in  the  GSC/NLMS  algorithm,  when  there  is  strong 
correlation  between  two  signals  arrivals  received  by  the  line  array. 

The  narrowband  outputs  of  the  GSC-SA/NLMS  algorithm  showed  that  in  this  case  the  signal 
cancellation  effects  have  been  minimized  and  the  two  signals  were  detected  only  at  the 
expected  two  beams  with  complete  cancellation  of  the  sidelobe  structure.  For  the  STMV 
beamformer,  the  grams  indicated  a  strong  sidelobe  structure  in  many  other  beams.  However, 
the  STMV-SA  beamformer  successfully  suppresses  the  sidelobe  structure  that  was  present  in 
the  case  of  the  STMV  beamformer.  From  all  these  simulations  [37],  it  was  obvious  that  the 
STMV-SA  beamformer,  as  a  narrowband  beamformer,  is  not  as  robust  for  narrowband 
applications  as  the  GSC-SA/NLMS. 
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Figure  11.  Concept  of  adaptive  sub-aperture  structure  for  line  arrays.  Schematic  diagram  shows  the 
steps:  (1)  formation  of  J  sub-apertures,  (2)  for  each  sub-aperture  formation  of  S  onventional  beams, 
(3)  fora  given  beam  direction,  6  formation  of  line  arrays  that  consist  of  J  number  of  directional 

sensors  (beams). 

5.2  Generic  Multi-Dimensional  Sub-Aperture 
Structure  for  Adaptive  Schemes 

The  decomposition  of  the  2-D  and  3-D  beamformer  into  sets  of  line  and/or  circular  array 
beamformers,  which  has  been  discussed  in  Section  3.2,  provides  a  first-stage  reduction  of 
the  numbers  of  degrees  of  freedom  for  an  adaptation  process.  Furthermore,  the  sub-aperture 
configuration  is  considered  in  this  study  as  a  second  stage  reduction  of  the  number  of 
degrees  of  freedom  for  an  adaptive  beamformer.  Then,  the  implementation  effort  for 
adaptive  schemes  in  multi-dimensional  arrays  is  reduced  to  implementing  adaptive  processes 
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in  line  and  circular  arrays.  Thus,  a  multi-dimensional  adaptive  beamformer  can  now  be 
divided  into  two  coherent  modular  steps  which  lead  to  efficient  system-oriented 
implementations. 


5.2.1  Sub-Aperture  Configuration  for  Line  Arrays 

For  a  line  array,  a  sub-aperture  configuration  includes  a  large  percentage 
overlap  between  contiguous  sub-apertures.  More  specifically,  a  line  array  is 
divided  into  a  number  of  sub-arrays  that  overlap,  as  shown  in  Figure  1 1 . 
These  sub-arrays  are  beamformed  using  the  conventional  approach;  and  this 
is  the  first  stage  of  beamforming. 

Then,  we  form  a  number  of  sets  of  beams  with  each  set  consisting  of  beams 
that  are  steered  at  the  same  direction  but  each  one  of  them  generated  by  a 
different  sub-array.  A  set  of  beams  of  this  kind  is  equivalent  to  a  line  array 
that  consists  of  directional  sensors  steered  at  the  same  direction,  with  sensor 
spacing  equal  to  the  space  separation  between  two  contiguous  sub-arrays  and 
with  the  number  of  sensors  equal  to  the  number  of  sub-arrays.  The  second 
stage  of  beamforming  implements  an  adaptive  scheme  on  the  above  set  of 
beams,  as  illustrated  in  Figure  1 1. 


Figure  12.  Concept  of  adaptive  sub-aperture  structure  for  circular  arrays,  which  is  similar  to  that  for 

line  arrays  shown  in  Figure  1 1. 

5.2.2  Sub-Aperture  Configuration  for  Circular  Array 

Consider  a  circular  array  with  M -sensors  as  shown  in  Figure  12.  The  first 
circular  sub-aperture  consists  of  the  first  M-G+l  sensors  with  n~l,2,...,M- 
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G+l,  where  n  is  the  sensor  index  and  G  is  the  number  of  sub-apertures.  The 
second  circular  sub-aperture  array  consists  of  M-G+l  sensors  with 
n=2,3,...,M-G+2.  The  sub-aperture  formation  goes  on  until  the  last  sub¬ 
aperture  consists  of  M-G+l  sensors  with  n=G,G+l,...,M.  In  the  first  stage, 
each  circular  sub-aperture  is  beamformed  as  discussed  in  Section  3.1.2  and 
this  first  stage  of  beamforming  generates  G  sets  of  beams. 

As  in  the  previous  section,  we  form  a  number  of  sets  of  beams  with  each  set 
consisting  of  beams  that  are  steered  at  the  same  direction  but  each  one  of 
them  generated  by  a  different  sub-array.  For  G<5,  a  set  of  beams  of  this  kind 
can  be  treated  approximately  as  a  line  array  that  consists  of  directional 
sensors  steered  at  the  same  direction,  with  sensor  spacing  equal  to  the  space 
separation  between  two  contiguous  sub-arrays  and  with  the  number  of  sensors 
equal  to  the  number  of  sub-arrays.  The  second  stage  of  heamforming 
implements  an  adaptive  scheme  on  the  set  of  beams,  as  illustrated  in  Figure 
12,  for  G=3. 

5.2.3  Sub-Aperture  Configuration  for  Cylindrical  Array 

Consider  the  cylindrical  array  shown  in  Figure  13  with  the  number  of  sensors 
X  =  NM ,  where  N  is  the  number  of  circular  rings  and  Mis  the  number  of 
sensors  on  each  ring.  Let  n  be  the  ring  index,  m  be  the  sensor  index  for  each 
ring  and  G  be  the  number  of  sub-apertures.  The  formation  of  sub-apertures  is 
as  follows: 

The  first  sub-aperture  consists  of  the  first  (N  -  G  + 1)  rings, 
where  n  =  1,2,...,N-G+1.  In  each  ring  we  select  the  first  set 
of  ( M  -G  +  l)  sensors,  where  m=l,2, ...M-G+l.  However, 
each  ring  has  Msensors,  but  only  (M-G+l)  sensors  are  used 
to  form  the  sub-aperture.  These  sensors  form  a  cylindrical 
array  cell,  as  shown  in  the  upper  right  hand  side  comer  of 
Fig.  13. 

In  other  words,  the  sub-aperture  includes  the  sensors  of  the  full  cylindrical 
array  except  for  G-l  sensors  from  G-l  rings,  which  are  denoted  by  small 
circles  in  Figure  13,  that  have  been  excluded  in  order  to  form  the  sub¬ 
aperture.  Next,  the  generic  decomposition  concept  of  the  conventional 
cylindrical  array  beamformer,  presented  in  Section  3.2.1,  is  applied  to  the 
above  sub-aperture  cylindrical  array  cell.  For  a  given  pair  of  azimuth  and 
elevation  steering  angles  { Gs ,  (j)s },  the  output  of  the  generic  conventional 
multi-dimensional  sub-aperture  beamformer  provides  beam  time  series, 
b  |  (t. ,  6s ,  (f)s ) ,  where  the  subscript  g=l  is  the  sub-aperture  index. 

The  second  sub-aperture  consists  of  the  next  set  of 

(N  —  G  + 1)  rings,  where  n  —  2,...,N-G+2.  In  each  ring  we 

select  the  next  set  of  (M  -G  +  l)  sensors,  where  m=2,...M- 
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G+2.  However,  each  ring  has  M  sensors,  but  only  (M-G+l) 
sensors  are  used  to  form  the  sub-aperture.  These  sensors  form 
the  second  sub-aperture  cylindrical  array  cell. 

Again,  the  generic  decomposition  concept  of  the  conventional  cylindrical 
array  beamformer,  presented  in  Section  3.2.1 ,  is  applied  to  the  above  sub¬ 
aperture  cylindrical  array  cell.  For  a  given  pair  of  azimuth  and  elevation 
steering  angles  { 6S ,  <ps },  the  output  of  the  generic  conventional  multi¬ 
dimensional  sub-aperture  beamformer  provides  beam  time  series, 
b  2  (fi  i  %  ’  0,  )  with  sub-aperture  index  g=2. 


The  proposed  sub-aperture  formation  continues  until  the  last 
sub-aperture  which  consists  of  a  set  of  (N  -  G  +  1)  rings, 
where  n  = G ,  G+1,...,N.  In  each  ring  we  select  the  last  set  of 
( M  —  G  + 1)  sensors,  where  m=G,G+l  ,...M.  Please  note  also 
that  each  ring  has  M  sensors,  but  only  (M-G+l)  sensors  are 
used  to  form  the  sub-aperture. 

As  before,  the  generic  decomposition  concept  of  the  conventional  cylindrical 
array  beamformer  is  applied  to  the  last  sub-aperture  cylindrical  array  cell.  For 
a  given  pair  of  azimuth  and  elevation  steering  angles  { 6 1  ,</>s },  the  output  of 
the  generic  conventional  multi-dimensional  sub-aperture  beamformer  would 
provide  beam  time  series,  b^=G(t  n  0s  ,  (j)s  )  with  sub-aperture  index  g=G. 

As  in  Section  5.2.2,  we  form  a  number  of  sets  of  beams  with  each  set 
consisting  of  beams  that  arc  steered  at  the  same  direction  but  each  one  of 
them  generated  by  a  different  sub-aperture  cylindrical  array  cell.  For  G<5,  a 
set  of  beams  of  this  kind  can  be  treated  approximately  as  a  line  array  that 
consists  of  directional  sensors  steered  at  the  same  direction,  with  sensor 
spacing  equal  to  the  space  separation  between  two  contiguous  sub-aperture 
cylindrical  array  cells  and  with  the  number  of  sensors  equal  to  the  number  of 
sub-arrays.  Then,  the  second  stage  of  beamforming  implements  an  adaptive 
scheme  on  the  above  set  of  beams,  as  illustrated  in  Figure  13. 

For  the  particular  case,  shown  in  Figure  13,  the  second  stage  of  beamforminu 
implements  an  adaptive  beamformer  on  a  line  array  that  consists  of  the  G-3 

beam  time  series  b^  (tnos,<pfi ) ,  g=l  ,2,...,G.  Thus,  fora  given  pair  of 

azimuth  and  elevation  steering  angles  {  0 v ,  (j)s },  the  cylindrical  adaptive 
beamforming  process  is  reduced  to  an  adaptive  line  array  beamformer  that 
includes  as  input  only  three  beam  time  series  b  (7(.,at,0v),  g-1,2,3  with 

spacing  8  —  [(R27T/M)2  +S1:]  ;  ,  which  is  the  spacing  between  two 
contiguous  sub-aperture  cylindrical  cells,  where  {R2jZ !  M )  is  the  sensor 
spacing  in  each  ring  and  S-  is  the  distance  between  each  ring  along  z-axis  of 
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the  cylindrical  array.  The  output  of  the  adaptive  beamformer  provides  one  or 
more  adaptive  beam  time  series  with  steering  centered  on  the  pair  of  azimuth 
and  elevation  steering  angles  {0S,  <j)s } . 

As  expected,  the  adaptation  process  in  this  case  will  have  near-instantaneous 
convergence  because  of  the  very  small  number  of  degrees  of  freedom. 
Furthermore,  because  of  the  generic  characteristics,  the  proposed  3-D  sub¬ 
aperture  adaptive  beamforming  concept  may  include  a  wide  variety  of 
adaptive  techniques  such  as  MVDR,  GSC  and  STMV  that  have  been 
discussed  in  References  [1,37], 


5.2.4  Sub-Aperture  Configuration  for  Planar  and  Spherical 
Arrays 

The  sub-aperture  adaptive  beamforming  concepts  for  planar  and  spherical 
arrays  are  very  similar  to  that  of  the  cylindrical  array.  In  particular,  for  planar 
arrays,  the  formation  of  sub-apertures  is  based  on  the  sub-aperture  concept  of 
line  arrays  that  has  been  discussed  in  Section  5.2.1.  The  different  steps  of 
sub-aperture  formation  for  planar  arrays  as  well  as  the  implementation  of 

adaptive  schemes  on  the  G  beam  time  series  bg  ( tj,6s ,  <j)s ) ,  g=l,2,...,G ,  that 

are  provided  by  the  G  sub-apertures  of  the  planar  array,  are  similar  with  those 
in  Figure  12  by  considering  the  composition  process  for  planar  arrays  shown 
in  Figure  7.  Similarly,  the  sub-aperture  adaptive  concept  for  spherical  arrays 
is  based  on  the  sub-aperture  concept  of  circular  arrays,  that  has  been 
discussed  in  Section  5.2.2. 
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FORMATION  OF  SUB-APERTURES  FORA 
MULTI-DIMENSIONAL  ARRAY  OF  SENSORS 

THE  SUB-APERTURE  FORMATION  IS  BASED  ON  DECOMPOSING  THE 
MULTI-DIMENSIONAL  BEAMFORMER  INTO  COHERENT  SETS  OF  LINE 
AND/OR  CIRCULAR  BEAMFORMERS 


GENERIC  STRUCTURE  DECOMPOSING  THE  3-D  BEAMFORMER 
OF  3-D  SUB-APERTURES  INTO  COHERENT  SUB-SETS  OF  LINE  8| 
CIRCULAR  ARRAY  BEAMFORMERS 

CONVENTIONAL  BEAMFORMING 

OUTPUT  PROVIDES  CONTINUOUS  BEAM  TIME  SERIES 


M1. 


A 


g=1,2 . G  IS  THE  SUB-APERTURE  INDEX 

i=1 ,2,..., I  IS  THE  INDEX  FOR  THE  TIME  SAMPLES| 
0  and  ARE  THE  AZIMUTH  AND  ELEVATION 
ANGLES  OF  STEERING 


9=1 


g=2 


g=3 


ADAPTIVE  BEAMFORMER 

FROM  THE  G  SUB-APERTURES  THE  BEAM  TIME  SERIES 
HAVING  THE  SAME  STEERING  DIRECTION  /  6  \ 

REPRESENT  G  DIRECTIONAL  SENSOR 
TIME  SERIES  OF  A  LINE  ARRAY 

PROVIDE  THE  ABOVE  DIRECTIONAL  TIME  SERIES  AT  THE  INPUT  OF  THE  ADAPTIVE  BEAMFORMER  AND 
GENERATE  S  NUMBER  OF  ADAPTIVE  BEAMS  WITH  STEERIN  DIRECTIONS  CENTERED  AT  (  0s><j>s)  for  s=1,2,...S 


OUTPUT  PROVIDES  CONTINUOUS 
ADAPTIVE  BEAM  TIME  SERIES 


Figure  14.  Signal  processing  of  generic  structure  decomposing  the  3-D  beamformer  for  cylindrical 
arrays  of  sensors  into  coherent  sub-sets  of  line  and  circular  array  beamformers. 
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5.3  Signal  Processing  Flow  of  a  3-D  Generic  Sub- 
Aperture  Structure 

As  it  was  stated  before,  the  discussion  in  this  report  has  been  devoted  in  designing  a  generic 
sub-aperture  beamforming  structure  that  will  decompose  the  computationally  intensive 
multi-dimensional  beamforming  process  into  coherent  sub-sets  of  line  and/or  circular  sub¬ 
aperture  array  bcamformers  for  ultrasound,  radar  and  integrated  active-passive  sonar 
systems.  In  a  sense,  the  proposed  generic  processing  structure  is  an  extension  of  a  previous 
effort  discussed  in  Reference  [1], 

The  previous  study  [1]  included  the  design  of  a  generic  beamforming  structure  that  allows 
the  implementation  of  adaptive,  synthetic  aperture  and  high-resolution  temporal  and  spatial 
spectral  analysis  techniques  in  integrated  active-passive  line-array  SONARs.  Figure  1 1  in 
Reference  [1]  shows  the  configuration  of  the  signal  processing  flow  of  an  equivalent  generic 
structure  for  line  arrays  that  allows  the  implementation  of  Finite  Impulse  Response  (FIR) 
filters  ,  conventional,  adaptive  and  synthetic  aperture  beamformers  [1,40,41,42], 

Shown  in  Figure  14  is  the  proposed  configuration  of  the  signal  processing  flow  that  includes 
the  implementation  of  line  and  circular  array  beamformers  as  Finite  Impulse  Response  (FIR) 
filters  [40,41,42].  The  processing  flow  is  for  3-D  cylindrical  arrays.  The  reconfiguration  of 
the  different  processing  blocks  in  Figure  14  allow's  the  application  of  the  proposed 
configuration  to  a  variety  of  ultrasound  systems  with  planar,  cylindrical  or  spherical  arrays 
of  sensors. 

As  discussed  at  the  beginning  of  this  section,  the  output  of  the  beamforming  processing 
block  in  Figure  14  provides  continuous  beam  time  series  for  ultrasound  image 
reconstruction.  Then  the  beam  time  series  are  provided  at  the  input  of  a  matched  filter  for 
Doppler  estimation.  This  modular  structure  in  the  signal  processing  ilow  is  a  very  essential 
processing  arrangement  in  order  to  allow  for  the  integration  of  a  great  variety  of  processing 
schemes  such  as  the  ones  considered  in  this  report.  The  details  of  the  proposed  generic 
processing  flow,  as  shown  in  Figure  14,  are  very  briefly  the  following: 

The  first  processing  block  includes  the  formation  of  sub-apertures  as  discussed  in  Section 
5.2.  Then  ,  the  sensor  time  series  from  each  sub-aperture  are  beamformed  by  the  generic 
multi-dimensional  beamforming  structure  that  has  been  introduced  in  Section  3.  Thus,  for  a 
given  pair  of  azimuth  and  elevation  steering  angles  {  0S ,  <j)s }  ,  the  output  of  the  generic 
conventional  multi-dimensional  beamformer  would  provide  Gbeam  time  series, 

0  (t  i,0s,  (fts )  ,  g—1,2 . G  .  The  second  stage  of  beamforming  includes  the  implementation 

of  an  adaptive  beamformer  as  discussed  in  Section  4.2. 
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6.  Concept  Demonstration:  Simulations  & 
Experimental  Results 


Concept  demonstration  of  the  proposed  ultrasonic  imaging  technology  has  been  carried  out 
with  a  custom  made  ultrasound  system  that  has  been  acquired  by  DRDC  to  fulfil  the 
project’s  development  requirements.  The  system  configuration  is  depicted  in  Figure  15.  The 
main  components  include: 

•  A  Line  array  of  32-transducers  for  ultrasound  imaging  applications. 

•  A  general  purpose  Ultrasound  Imaging  system  for  medical  diagnostic 
applications  with  CPU  components: 

•  Beamformer,  Data  manager  and  Display  functionality 

•  A  custom-made  A/DC  (12-bit,  45MHz  sampling  frequency  per  channel), 

•  a  signal  conditioning  unit  and  a  CPU  computer  workstation  that  includes 
implementation  of  the  advanced  beamforming  structure  discussed  in  the 
previous  sections. 


Figure  15.  Special  purpose  ultrasound  imaging  system  to  assess  image  resolution  improvements  of 

the  adaptive  beamformers. 


Thus,  the  sensor  time  series  (RF  data)  provided  by  the  line  array  (probe)  of  the  general 
purpose  ultrasound  system  of  Figure  14,  were  processed  with  the  advanced  beamforming 
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structure  of  the  present  investigation;  and  for  comparison  the  resulting  reconstructed  images 
were  compared  with  the  corresponding  images  of  the  general  purpose  ultrasound  system. 
Figure  16  presents  the  actual  custom  made  ultrasound  system  that  is  depicted  schematically 
also  in  Figure  15. 

In  summary,  the  basic  steps  of  the  concept  demonstration  process  included  the  following: 

•  The  sensor  time  series  (RF  data)  that  were  provided  by  the  line  array  of  the 
ultrasound  32-transducers  were  processed  by  the  general  purpose  ultrasound 
system  providing  as  an  output  a  complete  image  of  the  phantom  under 
investigation. 

•  The  same  RF  data  were  processed  also  by  the  computer  workstation  (CPU) 
including  the  adaptive  beamformers  (shown  in  Figure  15). 


Figure  1 6.  Left  picture  shows  the  general  purpose  ultrasound  system.  Right  picture  shows  the  data 
acquisition  system  providing  digitization  of  the  RF  data  from  the  probe  of  the  system  at  the  left. 

6.1  Computing  Architecture 

The  implementation  of  this  investigation’s  non-conventional  processing  schemes  in 
ultrasound  systems  is  a  non-trivial  issue.  In  addition  to  the  selection  of  the  appropriate 
algorithms,  success  is  heavily  dependent  on  the  availability  of  suitable  computing 
architectures. 

Past  attempts  to  implement  matrix-based  signal  processing  methods,  such  as  adaptive 
beamformers  reported  in  this  report,  were  based  on  the  development  of  systolic  array 
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hardware  because  systolic  arrays  allow  large  amounts  of  parallel  computation  to  be 
performed  efficiently  since  communications  occur  locally.  None  of  these  ideas  are  new. 
Unfortunately  systolic  arrays  have  been  much  less  successful  in  practice  than  in  theory.  The 
fixed  size  problem  for  which  it  makes  sense  to  build  a  specific  array  is  rare.  Systolic  arrays 
big  enough  for  real  problems  cannot  fit  on  one  board,  much  less  one  chip,  and  interconnects 
have  problems.  A  2-D  systolic  array  implementation  will  be  even  more  difficult.  So,  any 
new  computing  architecture  development  should  provide  high  throughput  for  vector  as  well 
as  matrix  based  processing  schemes. 

A  fundamental  question,  however,  that  must  be  addressed  at  this  point  is  whether  it  is 
worthwhile  to  attempt  to  develop  a  system  architecture  that  can  compete  with  a 
multiprocessor  using  stock  microprocessors.  Although  recent  microprocessors  use  advanced 
architectures,  improvements  of  their  performance  includes  a  heavy  cost  in  design 
complexity,  which  grows  dramatically  with  the  number  of  instructions  that  can  be  executed 
concurrently.  Moreover,  the  recent  microprocessors,  that  claim  high  performance  for  peak 
MFLOP  rates,  have  their  net  throughput  usually  much  lower  and  their  memory  architectures 
are  targeted  towards  general  purpose  code. 

The  above  issues  establish  the  requirement  for  dedicated  architectures,  such  as  in  the  area  of 
operational  sonar  and  radar  systems.  Sonar  applications  are  computationally  intensive  [1, 

60],  and  they  require  high  throughput  on  large  data  sets.  However,  the  experience  gained 
from  sonar  computing  architecture  developments  that  were  based  on  RISC  processors  [1], 
suggested  that  a  cost  effective  approach  in  that  direction  is  to  develop  a  multi-CPU 
computing  architecture  that  will  be  based  on  the  rapidly  evolving  microprocessor  technology 
of  the  CPUs  of  the  power-PCs,  such  as  the  family  of  Pentium  CPUs. 

Other  advanced  computing  architectures  that  can  cover  the  throughput  requirements  of 
computationally  intensive  signal  processing  applications,  such  as  those  discussed  in  this 
manuscript,  have  been  developed  by  Mercury  Computer  Systems.  Based  on  the  experience 
of  the  authors  of  this  report,  the  suggestion  is  that  implementation  efforts  of  advanced  signal 
processing  concepts  should  be  directed  more  on  the  development  of  generic  signal 
processing  structures  as  in  Figure  14,  rather  than  in  developing  very  expensive  computing 
architectures.  Moreover,  the  signal  processing  flow  of  advanced  processing  schemes  that 
include  both  scalar  and  vector  operations  should  be  very  well  defined  in  order  to  address 
practical  implementation  issues.  When  the  signal  processing  flow  is  well  established,  such  as 
in  Figure  14,  then  distribution  of  the  signal  processing  flow  in  a  number  of  parallel  CPU’s  of 
the  computing  architecture,  will  be  straightforward. 

In  the  following  sections,  we  address  the  practical  system  implementation  issues  by 
describing  the  current  effort  of  developing  a  real  3D  Ultrasound  system  deploying  a  planar 
array  to  address  the  requirements  of  the  Canadian  Forces  in  non-invasive  portable  diagnostic 
devices  deployable  in  operational  fields.  The  same  system  development  of  the  real  3D 
ultrasound  system  may  be  integrated  with  the  final  system  results  of  the  European 
Commission  project  “Minimally  Invasive  Tumour  Therapy  3D  Ultrasound  Guided” 
(MITTUG). 
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6.1.1  Experimental  Platform 


To  address  the  computing  architecture  requirements  of  this  investigation  we 
used  a  cluster  of  16  2-way  Pentium  III  nodes  interconnected  with  a  Myrinet 
network.  Myrinet  is  a  low-latency,  high-bandwidth,  point-to-point  system 
area  network  (SAN),  used  widely  for  clusters  of  workstations  and  PCs.  By 
allowing  users  to  directly  access  the  network,  without  operating  system 
intervention,  Myrinet  and  other  SANs  dramatically  reduce  latencies 
compared  to  traditional  TCP/IP  based  local  area  networks.  Moreover,  to 
further  reduce  latencies  in  SANs,  direct  memory  operations  are  usually 
supported;  reads  and  writes  to  remote  memory  are  performed  without  remote 
processor  intervention.  The  exact  system  configuration  is  summarised  in 
Table  1  and  for  more  details  the  reader  is  referred  to  [74],  The  nodes  in  the 
system  are  connected  with  a  Myrinet  interconnection  network.  Each  network 
interface  in  our  system  has  a  33  MHz  programmable  processor  and  connects 
the  node  to  the  network  with  two  unidirectional  links  of  160  MByte/s  peak 
bandwidth  each.  Actual  node-to-network  bandwidth  is  usually  constrained  by 
the  133  MBytes/s  I/O  bus  on  which  the  NIC  sits.  All  system  nodes  are 
connected  with  a  16-port  full  crossbar  Myrinet  switch. 

The  communication  layer  we  use  is  the  Message  Passing  Interface  (MPI)  on 
top  of  the  SCore  system  [75].  SCore  is  a  high-performance  parallel 
programming  environment  for  workstation  and  PC  clusters.  SCore  relies  on 
the  PMv2  [76]  low-level  communication  layer.  The  MPI  implementation  we 
use  is  a  port  of  the  MPICH  library  for  the  Score  system.  For  more  details 
about  the  bandwidth  and  latency  of  the  basic,  un-contended  MPI  Send  and 
MPI_Recv  operations,  the  reader  is  referred  to  [74], 


Table  1.  Cluster  node  configuration 


PROCESSORS 

2  x  Intel  Pentium  III,  800  MHz 

Cache 

Column  2  subheading 

Memory 

32K  (LI),  512K  (L2) 

OS 

512MB  SDRAM 

PCI  buses 

RedHat  Linux  Kernel  2.2.16-3smp 

NIC 

Myricom  M3M-PCI64M 

Communication  library 

MPICH/Score  4.0 

6.1.2  Algorithm  Implementation 

Implementations  of  the  signal  processing  structures  outlined  in  Sections  3  to 
5  assumes  that  input  data  is  delivered  directly  from  the  acquisition  unit  to  the 
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memories  of  the  nodes.  This  is  a  realistic  assumption,  since,  in  the  actual 
prototype  the  building  of  each  node  will  be  equipped  with  a  PCI-based  data 
acquisition  card  that  will  read  part  of  the  data  from  the  sensor  array.  Next,  we 
present  the  sequential  and  parallel  implementations  of  the  3D  beamforming 
algorithm. 

6.1 .2.1  Sequential  Implementation 

Our  sequential  algorithm  for  performing  the  computation  outlined 
in  Sections  3  to  5  and  outlined  in  Figure  14,  consists  of  the 
following  phases:  read  input  samples,  compute  FFTs,  filter  results, 
perform  column  steering,  reorganize  data  in  memory,  perform  row 
steering,  perform  inverse  FFTs,  and,  finally,  output  to  display. 

In  addition  to  dividing  the  scanned  volume  to  multiple  focal  zones 
and  using  multiple  beams  to  scan  it,  beamforming  algorithms 
divide  the  2D  surface  to  be  scanned  in  multiple  tiles.  If  we  view 
the  2D  surface  as  an  array  of  points,  we  can  divide  the  rows  and 
columns  in  ROW,  COL  groups  forming  ROW  x  COL  tiles.  Each 
tile  is  scanned  by  the  full  planar  sensor  array.  Thus,  channels-per- 
column  (CHC)  and  channels-per-row  (CHR)  represent  the  number 
of  sensors  in  each  dimension  of  the  planar  array.  Every  full 
snapshot  (that  generates  a  single  full  3D  frame  of  the  scanned 
volume)  requires  scanning  all  tiles  and  focal  zones,  as  illustrated 
in  Figure  17.  For  real-time  processing  we  would  require  at  least 
10  frames  (full  snapshots)  per  second  and  ideally  20  to  30. 

To  re-create  the  depth  information  the  algorithm  processes  the 
data  based  on  a  number  of  focal  zones  (NZONES).  Each  focal 
zone  is  of  constant  width  (depth)  which  depends  on  the  depth  of 
the  volume  to  be  scanned  and  the  number  of  focal  zones.  The 
volume  depth  is  usually  constant  and  defined  by  the  type  of 
objects  the  ultrasound  will  scan.  For  instance,  different  human 
organs  require  different  scan  depths.  In  this  work  we  use  a  fixed 
maximum  focal  depth  of  16  cm  and  we  vary  only  the  number  of 
focal  zones. 

For  each  scanned  point  of  the  input  volume  the  program  reads  the 
time  series  data  from  the  corresponding  sensor  (that  may  scan 
multiple  points)  and  stores  it  in  a  buffer  in  host  memory.  Each 
sample  is  32  bits  and  is  represented  as  a  single  precision  floating 
point  number. 

The  number  of  samples  that  need  to  be  processed  is  dictated  by 
the  depth  of  each  focal  zone.  The  probing  signal  used  to  scan  the 
object  reaches  different  depths  of  the  scanned  volume  with 
different  delays.  The  sampling  rate  used  to  digitize  the  received 
signal  dictates  the  minimum  number  of  samples  (and  the 
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minimum  FFT  size)  required  to  reconstruct  depth  information.  For 
example,  assuming  a  sampling  frequency  of  Fs  =  30  MFIz,  the 
number  of  samples  needed  to  reconstruct  20  cm  depth  can  be 
computed  by, 

N  =  Ml  =  2  x  0.2x30x10 6  =  7?92 
c  1540 

where  A  is  the  number  of  read  in  samples,  c  is  the  speed  of  sound 
in  meters  per  second,  d  is  the  depth  of  the  reconstruct  area  in 
meters.  The  factor  of  two  is  needed  to  account  for  the  round  trip 
time.  Thus,  in  this  example  the  ultrasound  system  (acquisition 
unit)  needs  to  provide  the  beamformer  with  7792  samples  for  each 
of  the  sensor  time  series  of  the  ultrasound  probe.  Using  more  than 
the  minimum  number  of  samples  can  improve  the  array  gain  and 
result  in  better  quality  imaging.  However,  reading  in  more 
samples  results  not  only  in  more  processing  but  also  in  longer 
acquisition  times  and  higher  storage  requirements.  To  counter  this 
problem,  in  our  experiments  we  set  the  number  of  acquired 
samples  to  8K  and  during  the  processing  we  vary  the  size  of  the 
FFT  operations. 

After  the  time  samples  are  read  and  converted  to  frequency 
domain,  a  filtering  phase  is  used  to  reduce  the  amount  of 
information  passed  to  later  stages.  The  information  embedded  in 
the  received  signals  necessary  for  reconstructing  a  particular  depth 
region  is  localized  in  a  certain  bandwidth.  Using  only  the  relevant 
frequency  components  further  reduces  computational  time.  Thus, 
the  FFT  output  samples  are  filtered  (with  a  Finite  Impulse 
Response  (FIR)  filter  [1])  to  exclude  unnecessary  information  and 
the  related  processing.  The  bandwidth  depends  on  the  focal  depth 
and  the  center  frequency  of  the  ultrasound  pulses.  Lower 
frequency  signals  usually  have  better  penetration  into  deeper 
regions,  whereas,  higher  frequency  signals  produce  sharper  beam 
resolutions.  However,  center  frequencies  are  usually  fixed  for 
each  depth.  Thus,  in  this  work  we  use  2  MHz  as  the  center 
frequency  for  an  input  volume  with  maximum  depth  of  16  cm.  In 
the  applications  we  are  interested  in,  most  objects  (human  organs) 
would  fall  within  this  range.  Given  this  center  frequency  the 
bandwidth  of  the  filter  can  vary  in  the  range  0. 5-4.0  MHz. 

After  filtering,  the  beamformer  performs  the  steering  operations 
and  finally  samples  are  converted  back  to  the  time  domain  for 
displaying.  Based  on  the  beam-steering  equations  (37),  the 
steering  of  beams  along  azimuth  and  elevation,  shown  in  Figure 
17,  separately  using  the  pre-calculated  steering  vector  to  align  the 
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time  delay  of  the  signals  arriving  in  different  sensors  and  IFFT 
transforms  the  signal  from  the  frequency  to  the  time  domain. 

6. 1.2.2  Optimiza tion 

To  gain  confidence  that  we  start  from  an  efficient  sequential 
implementation,  before  proceeding  with  parallelization,  we 
perform  a  number  of  measurements  to  fine  tune  several  aspects  of 
our  sequential  implementation. 

First,  we  explore  various  FFT  implementations,  both  our  own  and 
those  publicly  available.  We  find  that,  for  the  processors  we  use, 
the  most  efficient  implementation  is  FFTW  that  has  been 
developed  at  MIT,  a  C  library  for  computing  discrete  Fourier 
transforms,  which  is  free  for  non-commercial  use  under  the  GNU 
General  Public  License  (GPL)  [78],  This  also  minimizes  space 
requirements  since  the  output  of  this  function  is  a  half-complex 
array  that  consists  of  only  half  the  DFT  amplitudes;  The  negative- 
frequency  amplitudes  for  real  data  are  the  complex  conjugates  of 
the  positive-frequency  amplitudes.  The  side  effect  of  this  is  that 
we  need  to  re-organize  the  output  to  a  common,  full-complex 
array  format  after  the  FFT  and  revert  to  the  half-complex  array 
before  the  IFFT.  Also,  FFTW  computes  an  un-normalized 
transform  for  the  input  signal  IFFT(FFT(x))  =  N  x  X  for  size  N 
transforms.  Thus,  a  division  by  N  is  needed  for  each  element  of 
the  array  after  the  final  IFFT. 

Second,  we  experiment  with  multiple  ways  of  performing  the 
steering  and  the  related  dot  product  operations.  We  find  that  the 
best  results  are  obtained  by  using  the  function  from  the  Intel  Math 
Kernel  Library  (MKL)  [78]  to  compute  the  necessary  dot  products 
and  to  perform  the  steering.  MKL  is  optimized  for  the  Pentium 
family  of  processors  and  makes  effective  use  of  the  MMX  (Matrix 
Multiplication  Extensions),  SSE  (Streaming  SIMD  Extensions), 
and  similar  instructions. 

Third,  we  tune  loop  ordering  and  the  layout  of  multidimensional 
array  data  structures  to  improve  memory  access  behaviour  and  to 
reduce  cache  misses. 

The  overall  effect  of  these  optimization  steps  is  a  reduction  of  the 
overall  execution  time  of  the  sequential  implementation  by  a 
factor  of  about  10.  It  is  a  surprising  result  that  hand-tuning  can  be 
so  effective  with  all  compiler  optimizations  turned  on.  However, 
since  in  this  work  we  are  more  interested  in  the  behaviour  of  the 
parallel  version  we  omit  these  results  due  to  space  limitations. 


DRDC  Toronto  TR  2002-058 


67 


6. 1.2.3  Parallel  Implementation 

The  parallel  version  of  the  beamforming  algorithm  follows  closely 
the  structure  of  the  sequential  implementation.  We  see  that  the 
data  read  from  each  sensor  is  processed  independently  until 
steering.  Then,  during  the  steering  phase,  the  beams  across  the 
azimuth  and  elevation  directions  are  independent.  Therefore,  we 
choose  to  divide  the  computation  in  two  phases.  The  first  phase 
includes  all  processing  until  after  the  column-steering  phase.  The 
second  phase  includes  the  rest  of  the  processing,  starting  at  the 
row-steering  phase.  Between  the  two  phases,  we  need  to  re¬ 
organize  the  data  in  memory  by  performing  a  matrix  transpose, 
which  results  in  an  all-to-all  communication  pattern. 

The  first  phase  of  the  computation  for  each  frame  is  decomposed 
in  tasks  based  on  the  data  generated  by  each  sensor.  Thus,  there  is 
as  many  tasks  as  sensors  32x32  of  the  planar  array,  shown  in 
Figure  17,  which  is  sufficient  for  systems  with  large  numbers  of 
processors.  The  tasks  for  the  second  phase  are  determined  by  the 
processing  associated  with  each  beam.  We  use  the  processing 
related  to  a  single  beam  as  the  basic  task  and  we  decompose  the 
second  phase  to  (azimuth-BEAMS)x(elevation-BEAMS)  tasks. 
For  instance,  with  8  beams  in  each  direction,  there  are  64  coarse 
grain  tasks.  We  expect  that  for  all  practical  applications,  at  least 
8x8  beams  will  be  necessary  and  thus  we  do  not  consider  cases 
where  the  number  of  processors  is  larger  than  the  total  number  of 
beams. 

We  experiment  with  two  implementations  of  the  parallel 
algorithm.  The  first  implementation  uses  dedicated  nodes  for  each 
phase.  However,  we  find  that  balancing  the  number  of  nodes 
between  the  two  phases  of  the  computation  depends  on  a  large 
number  of  parameters.  Thus,  we  provide  a  second,  symmetric, 
implementation  as  well,  where  all  nodes  in  the  system  perform  the 
same  type  of  processing.  Although,  the  first,  dataflow  approach 
has  certain  advantages  in  reducing  task  management  costs,  we 
find  that  the  second  approach  is  more  flexible  and  results  in  better 
load  balancing.  Thus,  for  the  rest  of  this  work  we  only  use  our 
symmetric  implementation. 
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6.2  Details  of  a  Real  3D  Ultrasound  System  Deploying 
a  Planar  Array 

6.2.1  Requirements 

DRDC  Toronto’s  medical  imaging  project  investigations  include  the 
development  of  Ultrasound  Diagnostic  Imaging  Systems  with  adaptive  and 
synthetic  aperture  beamformers  for  improved  image  resolution  capabilities  to 
allow  for  tissue  identification.  The  system  will  be  capable  of  providing: 

2D  images  by  deploying  line  or  curvilinear  arrays, 

3D  volumetric  images  by  deploying  planar  arrays. 

The  main  R&D  activities  associated  with  this  investigation  include: 

1 .  Planar  Array  with 
32x32  -element  receiving  component 
12x1 2-element  transmitting  array 

The  corresponding  data  acquisition  system,  as  defined  in  Figure  18. 

2.  Development  of  a  PC-based  Computing  Architecture,  defined 
in  Figure  18,  to  provide  real-time  processing  of  raw  data 
provided  by  the  data  acquisition  system  defined  above.  The 
bandwidth  requirements  for  the  computing  architecture  are 
defined  in  Section  6.2.3  of  this  report. 

3.  Mapping  and  Implementation  of  the  synthetic  aperture  and 
adaptive  beamformers  in  the  dedicated  computing 
architecture  defined  above. 

4.  Development  of  the  synthetic  aperture  and  adaptive 
beamforming  structure  for  the  2D  and  3D  ultrasound  system 
of  this  investigation. 

5.  Development  of  3D,  4D,  6D  and  7D  image  visualization 
tools  for  the  output  beam  times  series  of  the  processing 
structure,  defined  in  the  above  activity  (4).  The  visualization 
tools  have  been  developed  by  Prof.  Sakas  at  Fraunhofer 
(National  Research  Institute,  Darmstadt,  Germany)  as  part  of 
our  group’s  involvement  in  the  European  Canadian 
collaborative  projects  New-Roentgen,  MITTUG,  ADUMS 
and  MRI-MARCB  [79-82], 
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6.  System  integration  of  R&D  efforts  of  the  above  activities, 
(1),  (2),  (3)  and  (4). 

In  the  following  sections,  the  system  characteristics  are  defined  in  terms  of 
the  acquisition  process,  bandwidth  of  data  flow  and  computing  architecture 
requirements. 


Figure  1 7.  Volume  4D  Digital  Scanning  for  Ultrasound  Applications 
using  Radar  &  Sonar  Phased  Array  Adaptive  Beamforming  of  DRDC 
Toronto  real  3D  ultrasound  system  development. 


6.2.2  Signal  Conditioning  and  Data  Acquisition  Unit 

The  signal  conditioning  unit  and  the  performance  characteristics  of  the  data 
acquisition  system  are  discussed  in  Section  6.2.3.  Briefly,  the  requirements 
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for  the  data  acquisition  process  of  the  received  ultrasound  energy  field  by  the 
omni-directional  receiving  transducers  of  the  32x32  planar  array,  include: 

The  digitization  process  of  the  32x8  sub-apertures  by  a  14-bit  256-channel 
A/D  unit  that  will  provide  the  signals  to  a  system  of  pin  connectors-cables 
with  suppressed  cross-talk  characteristics  (minimum  35  dB).  The  sampling 
frequency  is  33  MHz  for  each  of  the  channels  associated  with  a  receiving 
single  sensor.  The  multiplexer  associated  with  the  A/DC  allows  the  sampling 
of  the  32x32-sensors  of  the  planar  array  in  four  consecutive  active 
transmissions. 

The  computing  architecture,  (see  Table  1)  includes  sufficient  data  storage 
capabilities  for  the  sensor  time  series  as  indicated  by  the  system  requirements 
of  Figures  17  &  18.  The  A/D  and  signal  conditioning  modules  of  the  data 
acquisition  process  and  the  communication  interface  should  be  controlled 
with  software  drivers  that  form  an  integral  part  of  the  computing  architecture 
defined  in  the  following  section  6.2.3.  The  PC  based  computing  architecture 
includes  dedicated  software  that  allows  control  of  the  data  acquisition  and 
storage  process. 
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Figure  18.  Structure  of  the  PC-based  computing  architecture  development  for  the  real  3D  ultrasound  system  development. 


6.2.3 


Bandwidth  Requirements 


The  above  data  acquisition  characteristics  define  the  bandwidth  of  the  data  transfer 
process  as  follows: 

1.  The  volumetric  3D  ultrasound  system  will  cover  a  40°x40°  angular  sector  with 
maximum  penetration  depth  of  15  cm,  as  depicted  in  Figure  17.  The  transmitting 
beamformer  of  the  12x12  active  planar  array  steers  4  to  5  beams  across  azimuth  and 
4  to  5  beams  across  elevation.  This  process  defines  a  maximum  of  25  angular  sub¬ 
sectors,  each  of  10°xl0°  angular  size,  as  shown  in  Figure  17. 

2.  For  each  sub-sector  of  10°xl0°,  four  consecutive  active  transmissions  are  required  to 
illuminate  the  specific  sub-sector  of  interest  and  form  a  coherent  synthetic  aperture  of 
32x32  sensors. 

3.  For  a  penetration  depth  of  15  cm,  the  shortest  time  delay  between  two  consecutive 
active  transmissions  is  defined  by  (30/154000=)  x=0.2  ms. 

4.  Since  the  sampling  frequency  per  channel  of  the  A/D  unit  is  Fs  =  33  MHz,  the 
maximum  number  of  samples  for  each  receiving  time  series  is  8000  samples. 

5.  In  this  case,  the  requirement  for  data  transfer  is  (256-channels)  x  (4096-samples)  x 
(2-bytes)  =  2  Mb  in  0.2  ms 

6.  Thus,  for  a  specific  focus  range,  the  acquisition  process  requires: 

•  (4  transmissions)x(0.2  ms),  a  period  of  T=1.6  ms  to  sample  a  10°xl0°  sub-sector. 

•  In  this  case,  the  requirement  for  data  transfer  is  8x2  Mb  =  16  Mb  in  1.6  ms 

•  (9  transmissions)  x  (4  transmissions-per-sub-sector)  =  36  transmissions  or  a 
period  of  T=57.6  ms  to  sample  the  full  30°x30°  angular  sector. 

•  In  this  case,  the  requirement  for  data  transfer  is  9x16  Mb  =  144  Mb  in  57.6  ms. 

7.  To  track  cardiac  motion  effects  with  the  proposed  3D  US  system,  we  set  as 
requirement  to  achieve  10  samples  of  the  30°x30°  angular  sector  in  1  second. 
Although  the  number  of  10  samples  is  considered  small  to  track  most  of  the  heart’s 
complex  motion  effects,  this  limit  of  1 0  samples  for  4D  volumetric  visualization  of 
the  heart’s  motion  may  be  increased  at  a  later  time  without  major  modifications  in  the 
computing  architecture  structure  of  the  proposed  3D  ultrasound  system. 

8.  Thus,  the  requirement  for  data  transfer  for  4D  volumetric  imaging  is  10x144  Mb  = 
1.44  Gbytes  in  1  second  and  the  longest  time  delay  between  two  consecutive  active 
transmissions  may  be  as  long  as  t=1.587  ms,  to  accommodate  the  requirement  for  a 
total  number  of  360  transmissions  in  1  second. 

9.  If  the  computing  architecture  characteristics  of  the  data  acquisition  unit  will  allocate 
the  2  MB  of  data  from  each  one  of  the  4  transmissions  in  separate  memory  locations 
through  separate  bus-lines,  as  shown  in  Figure  18,  then  the  bandwidth  requirements 
are  defined  as  follows: 

10.  For  each  dedicated  bus-line,  the  requirement  of  data  transfer  is  reduced  to  1260/7 
Mbytes  =  180  Mbytes  in  1  second,  which  is  achievable. 

11.  Simultaneous  focusing  at  different  depths  is  achieved  by  designing  digital  pulses  that 
include  wavelets  with  different  frequency  regimes  and  different  frequency 
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modulation.  This  approach  is  conceptually  similar  with  that  being  used  in  sonar 
technology  [1], 

Figure  18  provides  a  schematic  overview  for  the  PC-based  computing  architecture  of 
DRDC  Toronto’s  3D  Ultrasound  system  under  development.  It  consists  of  two  major 
sections.  The  upper  section  defines  the  computing  architecture  for  the  data 
acquisition  unit.  The  lower  section  in  Figure  1 8,  provides  the  main  components  of  the 
PC-based  computing  architecture  for  the  processing  and  image  reconstruction  of  the 
ultrasound  sensor  time  series  that  has  been  briefly  reviewed  in  Section  6. 1 .  The  lower 
section  of  the  computing  architecture  is  divided  in  5  major  CPU  structures  that 
should  deliver  the  following  functionality: 

1)  CPU  for  synthetic  aperture.  The  7-processors  of  this  component,  shown  in 
Figure  18,  should  have  direct  connection  with  the  bus-lines  of  the  data 
acquisition  unit.  The  output  of  the  CPU  synthetic  aperture  structure  would 
provide  a  stream  of  (32x4096)x(32x4096)  time  series  for  adaptive 
beamforming.  Moreover,  the  memory  requirements  should  be  dedicated  to  this 
CPU  structure  as  local  memory. 

2)  CPU  for  the  transmitting  beamformer,  shown  at  the  left  hand  side  of  Figure  18. 
This  CPU  unit  estimates  the  active  steering  beams  (12xl2)x2048,  that  drives  the 
digital-to-analog  converter  (D/AC)  unit  of  the  active  sensors  of  the  planar  array. 

The  bandwhdth  requirements  for  the  transmitting  beamformer  are: 

•  9  sets  of  7  identical  beam  time  series,  each  of  size,  (12xl2)x2048x2  =  0.5  Mb 

•  Since,  for  a  given  focus  depth,  the  steered  beams  of  the  active  beamformer 
will  not  change  during  the  acquisition  process  of  the  10°xl0°  angular  sector 
in  1  second,  the  computing  architecture  of  the  data  acquisition  unit  should 
have  enough  memory  to  store:  9x0.5  Mb  =  4.5  Mb  of  data  for  a  given  focus 
range.  In  other  words,  for  a  given  range,  the  operation  of  the  3D  ultrasound 
system  will  initiate  a  data  transfer  of  4.5  Mb  to  the  memory  of  the  D/AC  unit 
of  the  data  acquisition  system,  before  the  actual  transmission  and  acquisition 
process  begins.  Thus,  the  bandwidth  requirements  in  this  case  are  very  small 
and  are  being  translated  into  memory  requirements  for  the  computing 
architecture  of  the  data  acquisition  unit,  associate  with  the  D/AC. 

3)  CPU  for  the  adaptive  beamformer,  which  is  shown  at  the  lower  part  in  Figure 
16.  This  CPU  sub-structure  requires  computationally  intensive  vector 
operations  of  size  (32x4096)xl  6  complex  numbers  =  2  Mb.  The  requirement  is 
for  several  vectors  to  be  processed  simultaneously,  which  leads  to  a  substantial 
parallelization  of  the  CPU  structure  with  significant  memory  requirements. 
DRDC  Toronto’s  adaptive  beamforming  structure  has  already  been  developed 
and  has  formed  the  basis  to  define  the  characteristics  of  the  CPU  structure  in 
terms  of  throughput  and  memory  requirements  discussed  in  Section  6.1. 

4)  The  output  of  the  CPU  adaptive  beamformer  provides  a  number  of  beams, 
generated  by  the  adaptive  beamforming  structure  with  the  following  beamwidth 
characteristics: 

•  It  has  been  assessed  that  the  proposed  adaptive  beamforming  structure 
provides  an  effective  beam-width  size,  which  is  equivalent  to  that  of  a  three 
times  longer  aperture  along  azimuth  and  three  times  longer  along  elevation  of 
the  deployed  planar  array.  Thus,  the  deployed  receiving  32x32  planar  array 


74 


DRDC  Toronto  TR  2002-058 


would  have  beamwidth  characteristics  equivalent  with  those  of  a 
(32x3)x(32x3)  size  planar  array. 

•  The  beam-width  of  a  receiving  32x32  planar  array  with  element  spacing  of 
0.5  mm  for  a  3  MHz  centre  frequency,  is  approximately  3.7°. 

•  Thus  the  receiving  adaptive  beams  along  azimuth  would  provide: 

•  (3  sub-angular-sectors)x(45x3)  =  405  beam  time  series  with  4096  samples 
each  along  azimuth  steering, 

•  (3  sub-angular-sectors)x(45x3)  =  405  beam  time  series  with  4096  samples 
each  along  elevation  steering. 

•  The  effective  angular  resolution  of  the  proposed  adaptive  system  in  terms 
of  beam-width  size  will  in  be  in  angular  sectors  of  less  than  1°  x  1°. 

•  The  above  angular  resolution  capabilities  of  the  adaptive  system  lead  to  the 
following  tissue  resolution  characteristics: 

•  For  C-scan  and  for  depth  of  10  cm,  the  1°  x  1°  angular  resolution  sector 
corresponds  to  a  (1.7  mm)  x  (1.7  mm)  =  2.89  mm2  size  of  tissue  resolution 
and  for  depth  of  5  cm  to  a  (0.8  mm)  x  (0.8  mm)  =  0.64  mm2  size  of  tissue 
resolution. 

•  For  B-scan  the  line  resolution  will  be  equivalent  to  the  wavelength  of  the 
transmitted  centre  frequency  ,  which  is  0.5  mm. 

•  Thus,  the  volume  resolution  of  the  3D  adaptive  beamforming  structure  will 
be  equivalent  to  (0.64  mm2 )  x  (0.5  mm)  -  0.32  mm3,  in  3D  tissue  size  at  a 
depth  of  5  cm. 

•  Finally,  the  output  of  the  adaptive  beamformer  will  provide  (405x405)x2048 
set  of  data  for  B-scan,  C-scan  and  S-scan  ultrasound  imaging  tomography 
applications  and  with  10  frames  per  second  to  track  dynamic  characteristics, 
especially  for  cardiac  imaging  applications. 

5)  The  CPU  structure  for  the  4D  volume  visualisation  of  the  output  of  the  beam¬ 
time  series  provides  the  post-processing  functionality  for  the  proposed  3D 
ultrasound  system. 

In  summary,  the  computing  architecture  of  the  proposed  system  has  the  following 
characteristics: 

Real  time  performance, 

PC  based, 

Parallelization  for  vector  based  signal  processing  including  digital  beamformers.  Figures 
16-17  and  Sections  5  and  6  provide  sufficient  details  to  assess  the  throughput 
requirements  associated  with  the  proj  ect  development  activity. 

The  A/DC  is  well  grounded  and  capable  to  sample  the  256  channels  with  an  equivalent 
14-bit  resolution  and  33  MHz  sampling  frequency  per  channel.  Moreover,  the  unit  has 
dedicated  memory  and  separate  bus-lines  to  achieve  the  data  transfer  rates  defined  above. 
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The  multiplexer  integrated  with  the  A/DC  allows  for  the  sampling  of  the  (32x32)  planar 
array  in  groups  of  seven  spatially  overlapping  sub-apertures  (8x32).  The  data  acquisition 
period  between  two  consecutive  sub-apertures  should  be  0.2ms. 

The  D/AC  is  capable  to  drive  144  channels  with  12-bit  resolution  and  33  MHz  sampling 
frequency.  The  period  between  two  consecutive  active  transmissions  is  in  the  range  of  0.2 
ms.  Moreover,  the  local  memory  of  the  D/AC  unit  has  the  capability  to  store  the  active 
beam  time  series  with  total  memory  size  of  4.5  Mb,  being  generated  by  the  main 
computing  architecture  for  each  focus  depth  and  transferred  to  the  local  D/AC  memory 
when  the  transmission-acquisition  process  begins. 
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6.3  Synthetic  Data 


6.3.1  Synthetic  Data  Results  for  Ultrasound  Systems  Deploying  Line 
Arrays 

Performance  assessment  and  testing  of  the  generic  conventional  and  sub-aperture 
adaptive  beamformers  that  have  been  discussed  in  this  report,  have  been  carried  out  with 
the  following  type  of  simulated  ultrasound  data: 

Synthetic  data  sets  that  have  been  generated  for  an  active  ultrasound  system  deploying  a 
line  array.  The  simulated  line  array  sensor  time  series  employed: 

•  An  aperture  size  of  linear  array:  N=32  sensors  with  detector  pitch  of 
0.4mm 

•  Sampling  Frequency  =  33  MHz 

•  Two  FM  pulses  of  duration  1 .52  (is. 

•  The  frequency  characteristics  defined  in  Table  2  for  the  first  three 
experiments. 

•  The  position  characteristics  defined  in  Table  2  for  the  first  three 
experiments. 

•  The  position  and  frequency  characteristics  defined  in  Table  3  for  the 
fourth  experiment. 


Table  2:  Parameters  for  Cases  I,  II  and  III. 


Point  Source  1  (Left) 

Point  Source  2  (Right) 

Fc 

BW 

Bearing 

Depth 

Fc 

BW 

Bearing 

Depth 

MHz 

MHz 

MHz 

MHz 

Case  I 

2.1 

1.1 

O 

oo 

65  mm 

2.0 

1.0 

96° 

65  mm 

Case  II 

3.6 

2.1 

0 

oo 

65  mm 

3.5 

2.0 

96° 

65  mm 

Case  III 

5.0 

3.0 

o 

oo 

65  mm 

5.0 

3.0 

96° 

65  mm 
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Table  3:  Parameters  for  Case  IV. 

Fc 

BW 

Bearing 

Depth 

Source  1 

4.0  MHz 

2.0  MFIz 

O 

o 

oc 

10  mm 

Source  2 

4.0  MHz 

2.0  MHz 

92° 

25  mm 

Source  3 

2.0  MHz 

1.0  MHz 

O 

T7f- 

OC 

40  mm 

Source  4 

2.0  MHz 

1.0  MHz 

96° 

50  mm 

The  synthetic  data  sets  generated  were  processed  with: 


•  Phased  array  conventional  beamformers  with  angular  sector  steering 
between  67.5°  and  1 1 2.5  °,  and  beam  resolution  of  0.57beam. 

•  Phased  array  sub-aperture  adaptive  STMV  beamformers  with  angular 
sector  steering  between  67.5°  and  1 12.5  °,  and  beam  resolution  of 
0.5°/beam. 


6. 3.1.1  Case  I:  Results 

The  resulting  broadband  beam  pattern  plots  from  the  experiments  with  the 
parameters  defined  as  Case  I  are  shown  in  Figure  19.  The  broadband  beam 
pattern  of  the  conventional  beamformer  with  spatial  shading  is  shown  by  the 
blue  curve,  and  from  the  conventional  beamformer  with  no  spatial  shading 
with  the  red  curve.  The  black  curve  shows  the  broadband  beam  pattern  for  the 
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STMV-SA  adaptive  beamformer.  These  plots  clearly  show  that  the  angular 
resolution  provided  by  the  STMV-SA  advanced  beamformer  is  superior  of  the 
conventional  beamformer.  The  green  and  orange  curves  show  the  broadband 
beam  patterns  derived  from  the  conventional  beamformer  with  spatial 
windowing  for  an  aperture  size  of  64  elements  and  96  elements  respectively. 

It  is  evident  that  the  effective  resolution  achieved  by  the  STMV-SA 
beamformer  with  32  elements  is  approximately  the  same  as  the  conventional 
beamformer  with  96  elements. 

Figure  20  shows  the  results  of  the  beamformers  in  the  time  domain.  The 
images  plot  the  envelopes  of  the  FM  signals  computed  using  the  Hilbert 
transform,  with  time  delay  (depth)  shown  vertically  and  beams  (angle)  shown 
horizontally.  The  top  row  of  images  from  left  to  right  show  the  reconstructed 
image  from  the  conventional  beamformer  with  32  detectors  and  spatial 
shading.  The  second  image  shows  the  result  from  the  64  element  conventional 
beamformer  without  spatial  shading.  The  third  image  shows  the  results  from 
the  conventional  beamformer  with  spatial  shading  for  a  96-element  array. 
These  images  correspond  to  the  blue,  red  and  orange  curves  in  Figure  19, 
respectively.  The  lower  row  shows  the  conventional  beamformer  with  spatial 
shading  for  a  32  element  on  the  left  and  the  STMV-SA  advanced  beamformer 
for  a  32-element  array  on  the  right.  These  images  correspond  to  the  green  and 
black  curves  of  Figure  19  respectively.  These  images  confirm  that  the  angular 
resolution  of  the  STMV-SA  beamformer  is  superior  and  equivalent  to  a  64 
element  conventional  beamformer  with  spatial  shading  and  approaches  the 
angular  resolution  of  a  96  element  conventional  beamformer  with  spatial 
shading. 
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Figure  20.  Reconstructed  Images  resulting 
from  Case  1 


Figure  21:  Reconstructed  Images  resulting 
from  Case  II. 


6. 3. 1.2  Case  II:  Results 

The  resulting  broadband  beam  patterns  from  the  experiments  with  the 
parameters  defined  as  Case  II  (see  Table  2),  are  shown  in  Figure  22.  The 
broadband  beam  pattern  of  the  conventional  beamformer  with  spatial  shading 
is  shown  by  the  blue  curve,  and  from  the  conventional  beamformer  with  no 
spatial  shading  with  the  red  curve.  The  black  curve  shows  the  broadband 
beam  pattern  for  the  STMV-SA  adaptive  beamformer.  These  plots  once  again 
show  that  the  angular  resolution  provided  by  the  STMV-SA  advanced 
beamformer  is  superior  of  the  conventional  beamformer  in  the  frequency 
regime  defined  for  Case  II.  The  green  and  orange  curves  show  the  broadband 
beam  patterns  derived  from  the  conventional  beamformer  with  spatial 
windowing  for  an  aperture  size  of  64  elements  and  96  elements  respectively. 
In  this  case,  the  effective  resolution  achieved  by  the  STMV-SA  beamformer 
with  32  elements  approaches  that  of  the  conventional  beamformer  with  64 
elements. 
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Beam  Cross  -  Section 


Figure  22.  Beam  pattern  resulting  from  Case  II. 

Figure  21  shows  the  results  of  the  beamformers  in  the  time  domain.  The 
images  plot  the  envelopes  of  the  FM  signals  computed  using  the  Hilbert 
transform,  in  the  same  manner  as  in  Figure  20.  The  top  row  of  images  from 
left  to  right  show  the  image  from  the  conventional  beamformer  with  32 
detectors  and  spatial  shading.  The  second  image  shows  the  result  from  the  64- 
element  conventional  beamformer  with  spatial  shading.  The  third  image 
shows  the  results  from  the  conventional  beamformer  with  spatial  shading  for  a 
96-element  array.  These  images  correspond  to  the  blue,  red  and  orange  curves 
in  Figure  19  respectively.  The  lower  row  shows  the  conventional  beamformer 
with  spatial  shading  for  a  32  element  and  the  STMV-SA  advanced 
beamformer  for  a  32-element  array.  These  images  correspond  to  the  green  and 
black  curves  of  Figure  19  respectively.  These  images  confirm  that  the  angular 
resolution  of  the  STMV-SA  beamformer  approaches  the  angular  resolution  of 
a  64  element  conventional  beamformer  with  spatial  shading. 

6. 3. 1.3  Case  III:  Results 

Figure  23  shows  the  resulting  broadband  beam  pattern  plots  from  the 
experiments  with  the  parameters  defined  as  Case  III  (see  Table  3).  The 
broadband  beam  pattern  from  the  conventional  beamformer  with  spatial 
shading  is  shown  by  the  blue  curve,  and  from  the  conventional  beamformer 
with  no  spatial  shading  with  the  red  curve.  The  black  curve  shows  the 
broadband  beam  pattern  for  the  STMV-SA  adaptive  beamformer.  Again,  in 
this  case,  the  angular  resolution  provided  by  the  STMV-SA  advanced 
beamformer  is  superior  to  the  conventional  beamformer  in  the  frequency 
regime  defined  by  the  parameters  of  Case  III  (see  Table  III).  However,  the 
improvement  is  less  marked  than  in  the  previous  two  cases.  The  green  and 
orange  curves  show  the  broadband  beam  patterns  derived  from  the 
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Figure  23.  Beam  pattern  resulting  from  Case  III. 
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Figure  24.  Images  resulting  from 

Figure  25.  Images  resulting  from  Case  IV. 

Case  III. 

Figure  24  shows  the  image  reconstruction  results  of  the  beamformers  in  the 
time  domain.  The  images  plot  the  envelopes  of  the  FM  signals  computed 
using  the  Hilbert  transform,  in  the  same  manner  as  in  Figures  20  and  21. 
Since  the  frequency  regime  of  Case  III  is  higher  than  the  design  frequency  of 
the  array,  this  affects  the  STMV-SA  beamformer’s  performance;  and  the 
gains  achieved  are  not  as  evident  as  in  the  previous  two  cases.  Therefore,  the 
conclusion  from  these  simulations  is  that  for  the  advanced  beamformers  to 
achieve  optimum  performance,  in  terms  of  angular  resolution,  it  is  essential 
that  the  frequency  regime  of  the  deployed  ultrasound  signals  should  be  close 
to  the  design  frequency  of  the  array. 


6. 3. 1.4  Case  IV:  Results 

The  images  in  Figure  25  show  the  operation  of  the  beamformers  in  a  more 
general  case.  The  parameters  of  the  sources  are  defined  in  Table  3.  The  upper 
two  images  show  the  results  from  the  conventional  beamformer  with  spatial 
window  for  a  32-element  array  on  the  left,  and  a  64-element  array  on  the 
right.  The  lower  two  images  show  the  conventional  beamformer  with  no 
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spatial  shading  for  a  32-element  array  on  the  left  and  a  64  element  array  on 
the  right.  The  rightmost  image  shows  the  results  from  the  STMV-SA 
beamformer  for  a  32-element  array.  Comparison  of  the  results  from  the 
STMV-SA  beamformer  and  the  conventional  beamformer  with  spatial 
shading  show  that  the  STMV-SA  beamformer  provides  better  angular 
resolution  than  the  conventional  beamformer  employing  a  64  element  array. 
Without  spatial  shading,  the  angular  resolution  of  the  64  element  conventional 
beamformer  is  approximately  equivalent  with  that  of  the  STMV-SA 
beamformer,  but  it  suffers  from  sidelobe  effects,  which  are  most  evident  with 
the  closest  source.  Overall,  these  results  show  that  the  STMV-SA  beamformer 
provides  improved  angular  resolution  over  the  conventional  beamformer  with 
spatial  shading.  In  addition,  it  provides  improved  sidelobe  suppression  over 
that  of  the  conventional  beamformer  with  no  spatial  shading. 


6.3.2  Synthetic  Data  Results  for  Ultrasound  Systems  Deploying  Planar 
Arrays 

The  synthetic  data  experiments  for  the  planar  array  was  carried  out  using  the  Field  II 
ultrasound  simulator  program  obtained  from  the  Technical  University  of  Denmark 
('http://www.it.dtu.dk/~iai/fieldl.  The  simulator  simulates  point  sources,  and  was  set  up  to 
simulate  5,000  point  sources  arranged  in  a  spherical  shell,  conforming  to  the 
specifications  of  the  system  defined  in  Section  6.2.  Specifically,  a  32x32  planar  array  on 
receive  and  a  12x12  array  on  transmit  were  simulated,  with  element  spacing  of  0.3  mm 
and  sampling  frequency  of  33  MHz.  The  FM  pulse  was  centred  at  2.5  MHz,  with  a 
bandwidth  of  4.0  MHz.  The  illumination  pattern  is  as  described  in  Figure  17,  but  with  six 
transmit  beams  spaced  10°  apart  covering  60°  is  each  direction  (azimuth  and  elevation). 
The  result  was  a  total  of  36  sectors.  A  conventional  beamformer  was  used  to  process  the 
data  received  by  the  32x32  array.  The  array  decomposition  was  preformed  as  defined  in 
Section  5.2.4,  processing  each  row  of  the  array  to  obtain  azimuth  beams,  and  then 
processing  the  azimuth  beams  for  a  given  direction  for  all  of  the  rows  to  obtain  3D 
azimuth  and  elevation  beams.  A  complete  reconstruction  of  the  beams  into  a  volume  was 
then  performed. 

The  C-scans  derived  from  the  3D  reconstructed  volumes  of  the  spherical  shell  are  shown 
in  Figure  26.  In  this  image,  which  shows  a  slice  of  the  spherical  shell,  the  expected  ring 
that  corresponds  to  the  cross  section  of  a  shell  is  visible.  The  left  and  right  images  in 
Figure  26  correspond  to  the  3D  conventional  and  adaptive  beamformers,  respectively. 

The  better  performance  of  the  adaptive  beamformer  is  evident  in  this  case  as  it  provides 
better  detection  and  image  definition  that  the  corresponding  conventional  beamformer. 
The  visualization  software  was  provided  by  Prof.  Sakas  of  Fraunhofer  (National  Research 
Institute,  Darmstadt,  Germany)  [60],  as  part  of  our  technical  exchanges  within  the 
framework  of  the  European-Canadian  project  MITTUG  [81]. 
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Figure  26.  C-scans  derived  from  the  3D  reconstructed  images  of  the  simulated  spherical  shell.  Left  image. 
reconstructed  with  the  3D  conventional  beamformer.  Right  image,  reconstructed  with  the  3D  adaptive 

beamformer. 

Figure  27  shows  the  3D  volume  reconstruction  of  the  spherical  shell  using  the  3D 
conventional  (left  image)  and  the  3D  adaptive  (right  image)  ultrasound  beamformers 
defined  in  the  present  report.  As  was  the  case  with  the  C-scans  (Figure  26),  the  results  in 
Figure  27  show  that  the  3D  adaptive  beamformer  provides  better  image  definition  than 
the  corresponding  3D  conventional  beamforming  results. 


Figure  27.  3D  volume  reconstruction  of  the  simulated  spherical  shell.  Left  image,  reconstructed  with  the  3D 
conventional  beamformer.  Right  image,  reconstructed  with  the  3D  adaptive  beamformer. 
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6.3.3  Real  Data  Results  for  Ultrasound  Systems  Deploying  Linear 
Arrays 

The  real  data  tests  were  carried  out  using  the  equipment  described  at  the  beginning  of  this 
section  (Section  6)  and  shown  in  Figures  15  and  16.  The  ultrasound  system  employs  a  96- 
element  array,  of  which  only  32  are  active  at  any  given  time.  The  first  32  elements  of  the 
array  are  selected,  and  8  of  these  elements  transmit  an  RF  pulse  to  illuminate  the  object 
of  interest.  The  illumination  is  focused  to  the  broadside  of  the  32  elements  selected,  and 
at  a  selected  focal  distance.  The  32  elements  then  receive  reflections.  This  is  for  one 
active  transmission.  For  the  next  transmission,  a  new  set  of  32  elements  is  selected  and 
the  process  is  repeated.  The  new  32  elements  are  selected  in  a  sliding  window  fashion, 
using  31  from  the  previous  active  transmission,  and  one  new  element.  In  other  words,  the 
active  elements  slide  over  (or  pitch)  by  one  detector.  This  process  in  repeated  by  sliding 
one  detector  for  each  transmission  until  64  active  transmissions  are  made,  and  all 
elements  have  been  exhausted.  The  experiments  were  carried  out  using  data  acquired  in 
this  manner.  In  the  standard  operation  of  the  ultrasound  system  more  active  illuminations 
are  achieved  by  using  the  so  called  “edge-effects”  where  the  active  elements  begin  with 
16  active  elements  initially,  and  increase  with  each  pitch  (or  active  transmission)  to  a 
maximum  of  32.  When  the  last  valid  pitch  that  allows  32  active  detectors  is  reached,  the 
system  continues  until  decreasing  the  number  of  active  elements  until  there  are  only  16 
elements  active.  In  other  words,  in  standard  operation,  the  system  pretends  that  there  are 
16  non-  working  detectors  before  the  96-element  array,  and  16  more  after  it,  making  it  a 
potential  128-element  array.  By  using  a  sliding  window  of  length  of  32  elements,  then 
there  are  96  active  transmissions  (and  reflections  received)  possible.  In  the  experiments, 
only  the  64  pitches  where  32  physical  elements  were  active  were  used,  hence  there  are  64 
transmissions  to  broadside.  The  received  reflections  are  beamformed  to  broadside,  and  a 
selected  focal  depth  for  these  64  data  segments  received.  This  yields  an  image  with  64 
rows  -  one  for  each  active  transmission  position. 

The  phantoms  used  in  the  experiments  were: 

Two  wire  cross  sections  simulating  two  point  sources  in  space  suspended  in  a  water  bath. 
A  ball  suspended  in  a  water  bath. 

The  resulting  images  from  the  ultrasound  system,  the  conventional  beamformer  with 
spatial  shading,  and  the  STMV-SA  beamformer  are  shown  in  Figures  28  and  29.  The 
leftmost  image  of  Figure  28  show's  the  image  obtained  from  the  ultrasound  system  under 
standard  operating  conditions.  The  centre  image  show's  the  results  obtained  by  processing 
the  raw  data  captured  by  the  data  acquisition  with  a  conventional  beamformer  including 
spatial  shading.  The  right  most  image  show's  the  image  obtained  by  processing  the  raw 
data  with  the  STMV-SA  beamformer.  These  two  images  to  the  right  are  plotted  using  a 
logarithmic  scale  to  show  the  background  effects  more  clearly.  The  left  image  and  the 
two  right  images  do  not  employ  the  same  magnification  factor.  Comparison  of  these 
images  shows  the  improvement  in  detection. 

There  are  three  major  factors  that  affect  the  performance  of  the  beamformers  in  this 
scenario.  First,  the  frequency  regime  is  W'ell  outside  of  the  design  frequency  of  the  array. 
With  an  element  spacing  of  approximately  0.4  mm  the  design  frequency  of  the  array  is 
approximately  1.875  MFIz,  while  the  centre  frequency  of  the  RF  pulses  is  approximately 
8  MHz  with  8  MHz  bandwidth.  The  second  factor  is  the  lack  of  dynamic  range  of  the 
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A/D  converter,  which  is  equivalent  to  4-bit.  Thus,  the  A/D  converter  does  not  exploit  the 
full  dynamic  range  during  data  acquisition,  thereby  compromising  the  SNR  of  the  system. 
The  third  factor  is  the  sensitivity  of  the  data  acquisition  system  to  background 
interferences.  The  interference  pattern  created  by  this  in  combination  with  the  lack  of 
dynamic  range  is  visible  as  diagonal  lines,  and  spots  that  form  a  diagonal  pattern  in 
Figures  28  and  29.  Figure  29  shows  the  results  with  the  ball  phantom  in  the  same  format 
as  Figure  28.  Despite  the  interference,  the  improvement  in  detection  by  the  STMV-SA 
beamformer  is  evident.  There  are  currently  efforts  underway  to  improve  the  cabling  and 
electronics  of  the  data  acquisition  system,  shown  in  Figure  16,  to  reduce  the  sensitivity  to 
interference.  In  addition,  variable  amplifiers  will  be  added  to  exploit  the  full  dynamic 
range  of  the  system’s  A/DC  peripheral. 


Figure  28.  Results  with  wire  phantom. 


DRDC  Toronto  TR  2002-058 


87 


Figure  29.  Results  with  Ball  Phantom. 
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7.  Conclusion 


The  synthetic  and  real  data  results  of  this  report  indicate  that  the  generic  multi-dimensional  adaptive 
concept  addresses  practical  concerns  of  near-instantaneous  convergence,  shown  in  Figures  18-29,  for 
ultrasound  systems.  The  performance  characteristics  of  the  sub-aperture  adaptive  beamformer  compared 
with  that  of  the  conventional  beamformer  are  reflected  as  improvements  in  directional  estimates  of 
azimuth  and  elevation  angles  and  suppression  of  reverberation  effects.  This  type  of  improvement  in 
azimuth  and  elevation  bearing  estimates  is  essential  for  3-D  ultrasound  operations. 

In  summary,  a  generic  beamforming  structure  has  been  developed  for  multi-dimensional  sensor  arrays 
that  allows  the  implementation  of  conventional,  synthetic  aperture  and  adaptive  signal  processing 
techniques  in  integrated  active-passive  real  time  systems.  The  proposed  implementation  is  based  on 
decomposing  the  2-D  and  3-D  beamforming  process  in  sub-sets  of  coherent  processes  and  creating  sub¬ 
aperture  configurations  that  allow  the  minimization  of  the  number  of  degrees  of  freedom  of  the  adaptive 
processing  schemes.  The  proposed  approach  has  been  applied  to  line,  planar,  and  cylindrical  arrays  of 
sensors  where  the  multi-dimensional  beamforming  is  decomposed  into  sets  of  line-array  and/or  circular 
array  beamformers.  Moreover,  the  application  of  spatial  shading  on  the  generic  multi-  dimensional 
beamformer  is  a  much  simpler  process  than  that  of  the  fully  coherent  3-D  beamformer.  This  is  because 
the  decomposition  process  allows  two  simple  and  separate  applications  of  spatial  shading  (i.e.  one  for 
circular  and  the  other  for  line  arrays). 

The  beamformers  have  been  applied  to  real  and  simulated  data  sets  of  ultrasound  systems.  In  general 
ultrasound  systems  require  near-field  beamforming.  The  results  with  the  real  and  synthetic  data  sets 
show  that  firstly,  the  decomposition  and  sub-aperture  techniques  hold  true  for  focused  beamformers. 
Also,  the  results  show  that  the  performance  of  the  advanced  beamformer  is  not  affected  by  the  near-field 
conditions  presented  in  ultrasound  systems,  and  that  performance  gains  over  conventional  beamformers 
are  achieved. 

The  fact  that  the  sub-aperture  adaptive  beamformers  provided  array  gain  improvements  for  FM  signals 
under  a  real  time  data  flow  as  compared  with  the  conventional  beamformer  demonstrates  the  merits  of 
these  advanced  processing  schemes  for  practical  ultrasound  applications.  In  addition,  the  generic 
implementation  scheme  of  this  study  suggests  that  the  design  approach  to  provide  synergism  between 
the  conventional  beamformer,  the  synthetic  aperture  and  the  adaptive  processing  schemes,  (e.g.,  see 
results  in  Figures  22  and  23)  is  an  essential  property  for  system  applications. 

Although  the  focus  of  the  implementation  effort  included  only  a  few  adaptive  processing  schemes,  the 
consideration  of  other  types  of  spatial  filters  for  real  time  ultrasound  applications  should  not  be 
excluded.  The  objective  here  was  to  demonstrate  that  adaptive  processing  schemes  could  address  some 
of  the  challenges  that  the  next  generation  ultrasound  systems  will  have  to  deal  with  in  the  near  future. 
Once  a  generic  signal  processing  structure  is  established,  as  suggested  in  the  report,  the  implementation 
of  a  wide  variety  of  processing  schemes  can  be  achieved  with  minimum  efforts  for  real  time  systems 
deploying  multi-dimensional  arrays.  Finally,  the  results  presented  in  this  report  indicate  that  the  sub¬ 
aperture  STMV  adaptive  scheme  address  the  practical  concerns  of  near-instantaneous  convergence 
associated  with  the  implementation  of  adaptive  beamformers  in  ultrasound  systems. 
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9.  List  of  symbols  /  abbreviations  /  acronyms  /  initialisms 


0* 

A/DC 

4(f<) 

a 

AG 

BW 

KfA) 


B(f,0s) 


B(0) 

CFAR 

C 

c 

DT 

5 

d 


£ 


Complex  Conjugate  Transpose  Operator 
Analog-to-Digital  Converter 
Power  spectral  density  of  signal  s(tj) 

Small  positive  number  designed  to  maintain  stability  in  NLMS  adaptive  algorithm 
Array  Gain 
Signal  Bandwidth 

Beams  for  conventional  or  adaptive  beamformers,  of  an  array  steered  at  direction  9S  and 
expressed  by  />  ) =  £  (/Mi  (/>  ^ ) 

Narrowband  beam  power  pattern  of  an  array  expressed  by 

B(f,ej  =  b(f,es)b'(f,es) 

Broadband  beam  power  pattern  of  an  array  steered  at  direction  9 

Constant  False  Alarm  Rate 

Signal  blocking  matrix  in  GSC  adaptive  algorithm, 

speed  of  sound  in  a  medium  (i.e.  human  body) 

Detection  Threshold 

sensor  spacing  for  a  line  array  receiver 

Detection  Index  of  Receiver  Operating  Characteristic  (ROC)  curve 

Steering  vector  having  its  nth  phase  term  for  the  signal  wave  arrival  with  angle  9  being 

expressed  by  dn  (Jt,  6)  =  expjy‘2  K'1  T„(0)  j  , 

Noise  vector  component  with  n  th  element  £n(tj)  for  sensor  outputs 

(i.e.  x  =  S  £) 

Expectation  operator 


DRDC  Toronto  TR  2002-058 


95 


ETAM 

e 

f 

fs 

®M,) 

GSC 

h 

1 

i 

k 

k 

LCMV 

1 

L 

MVDR 

M 

M 

N 

NLMS 

Ne 


Extended  Towed  Array  Measurements 

Angle  of  plane  wave  arrival  with  respect  to  an  array  receiver 

Frequency  in  Hz 

Sampling  frequency 

Steered  spatial  covariance  matrix  (STCM)  in  time  domain 
Generalized  Sidelobe  Canceller 

Vector  of  weights  for  spatial  shading  in  beamforming  process 
Unit  vector  of  ones 

Index  of  time  samples  of  sensor  time  series,  {x„(tj)  ,i=l,2,...,M} 

wavenumber  parameter 

Iteration  number  of  adaptation  process 

Linear  Constraint  Minimum  Variance 

wavelength  of  signal  with  frequency/  where  c-fX 

Size  of  line  array  expressed  by  L—  (N-l)S 

Minimum  Variance  Distortionless  Response 

Number  of  time  samples  in  sensor  time  series,  where  {x„(tj)  ,i=l,2,...,M } 

Convergence  controlling  parameter  or  “step  size”  for  the  NLMS  algorithm. 

Number  of  sensors  in  an  array  receiver,  where  {x„(tj)  ,n=l,2,...,N  } 

Normalized  Least  Mean  Square 

Noise  energy  flux  density  at  the  receiving  array 

Index  for  space  samples  of  sensor  time  series  {x„(tj)  ,n=],2,...,N } 
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N 

<^(0  >ti )  Beam  time  series,  output  of  time-domain  beamformer,  £(ds ,  tx )  =  ^  xn  (7;  -Ts)  or 

n=  1 

formed  by  using  FFTs  and  fast  convolution  of  frequency  domain  beamformers,  £(9s>tj)  = 
IFFT{^^} 

OMI  Operator-Machine  Interface 

P]3  Probability  of  Detection  for  ROC  curves 

PfA  Probability  of  False  Alarm  for  ROC  curves 

n  3.14159 

R(fi )  Spatial  correlation  matrix  with  elements  R„m(f,d„m)  for  received  sensor  time  series 

Pnm(f'8nm)  Crosscorrelation  coefficients  given  from,  p„m  (f  ’dnm)  =  Rnm(f,dnm)/X2(f) , 

ROC  Receiver  Operating  Characteristic  (ROC)  curve 

S  signal  vector  whose  nth  element  is  expressed  by  sn(ti)  =  sn[ti  +  T  n(Q)\ 

S  Spatial  correlation  matrix  for  the  received  signal  sn(tj) 

S(fi,  0)  Spatial  correlation  matrix  for  the  plane  wave  signal  in  Frequency  domain.  It  has  as  its  nth 

row  and  nzth  column  defined  by,  Snm(fv  0)  =  As(fi)dn(fi,  d)d*m(flt  6). 

STCM  Steered  Covariance  Matrix 

STMV  Steered  Minimum  Variance 

SVD  Singular  Value  Decomposition  method 

<J2n  (  f  )  )  Power  spectral  density  of  noise,  £n(lj). 

X  Row  vector  of  received  N-  sensor  time  series  {xn(tj)  ,n=l,2,...,N } 

Xn  (f)  is  the  F ourier  transform  of xn(tj) 

X  n  (f. ,  6S  )  Pre-steered  sensor  time  series  in  frequency  domain.  In  time  domain  are  denoted  by 
xn  {ti ,  Tn  (6s ))  Pre-steered  sensor  time  series  in  time  domain 
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X2{f ) 

T(fk,9) 


Mean  acoustic  intensity  of  time  sequences  at  frequency  bin/ 

Time  delay  between  the  1st  and  the  nth  sensor  line  array  for  an  incoming  plane  wave  with 
direction  of  propagation  0 , 

Diagonal  steering  matrix  with  elements  those  of  the  steering  vector,  D^f^d) 


Wifi, 9) 
co 


Adaptive  steering  vector 
Frequency  in  radians/second 

Result  of  the  signal  blocking  matrix  C  being  applied  to  presteered  sensor  time  series 


x(f„e,) 
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14.  ABSTRACT 

(U)  The  present  report  outlines  the  development  of  an  advanced  digital  ultrasound  imaging  technology 
leading  to  a  next-generation  field-portable  4D  diagnostic  imaging  system.  The  proposed  technology 
consists  of: 

„h  Adaptive  beamforming  schemes  [1-6]  to  provide  high  image-resolution  for  ultrasound  diagnostic 
systems 

„h  3D  visualization  imaging  techniques  to  assist  imaging  during  field-deployable  minimally  invasive 
operations 

„h  time-reversal  pulse  design  for  ultrasound  signals  to  eliminate  aberration  effects,  which  cause  fuzziness 
in  reconstructed  ultrasound  images;  and 

„h  investigations  on  computing  architectures  and  planar  ultrasound  sensor  arrays  allowing  system 
integration  of  the  proposed  technologies  as  a  compact-portable  system  for  field-deployable  operations. 
The  aim  is  to  investigate  a  new  ultrasonic  sensing  and  imaging  technology  that  would  lead  to  the  design 
of  portable,  compact  and  field  deployable  3D  ultrasound  diagnostic  systems.  The  performance 
characteristics  of  the  proposed  diagnostic  systems  include  minimization  of  the  aberration  effects  and 
significant  improvement  in  image  resolution  to  allow  for  tissue  identification  and  3D  tomographic 
imaging  of  internal  body-organs.  The  adaptive  processing  schemes  of  this  study  are  characterized  as 
Dual-Use  technologies  and  have  been  tested  in  operational  active  sonars  of  the  Canadian  Navy  [1,4-6]. 
Real  data  results  have  shown  that  they  provide  array  gain  improvements  for  signals  embedded  in 
anisotropic  noise  fields,  similar  to  those  in  the  human  body. 
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